Extract Shiny90 output and aggregate for cryptococcal meningitis CD4 <200 inputs

ahd-first90-sim_2024-04-27.R

library(dplyr)
library(brio)
library(readr)
library(first90)  # remotes::install_github("mrc-ide/first90release")


#' ## Define functions

get_pjnz_shiny90_filename <- function (pjnz) {
  files <- utils::unzip(pjnz, list = TRUE)$Name
  shiny90file <- grep("\\.shiny90$", files, ignore.case = TRUE, value = TRUE)
  if (length(shiny90file) > 1) {
    msg <- paste0("Multiple .shiny90 files found: ", paste0(shiny90file, collapse = ", "))
    shiny90file <- shiny90file[which.min(nchar(shiny90file))]
    msg <- paste0(msg, "\nUsing file: ", shiny90file)
    warning(msg)
  }
  shiny90file
}

get_proj_years <- function(ss) {
  ss$proj_start + 1:ss$PROJ_YEARS - 1L
}

first90_mod_add_dimnames <- function (mod, ss) {
  yrlbl <- get_proj_years(ss)
  sexlbl <- c("male", "female")
  hAGlbl <- c("15-16", "17-19", "20-24", "25-29", "30-34",
              "35-39", "40-44", "45-49", "50+")
  pAGlbl <- 15:80
  pDSlbl <- c("hivn", "hivp")
  hDSlbl <- c(">500", "350-499", "250-349", "200-249", "100-199",
              "50-99", "<50")
  hTSlbl <- c("art0mos", "art6mos", "art1yr")

  dimnames(mod) <- list(age = pAGlbl, sex = sexlbl, hiv = pDSlbl,  year = yrlbl)
  dimnames(attr(mod, "hivpop")) <- list(cd4stage = hDSlbl, agegr = hAGlbl, sex = sexlbl, year = yrlbl)
  dimnames(attr(mod, "diagnpop")) <- list(cd4stage = hDSlbl, agegr = hAGlbl, sex = sexlbl, year = yrlbl)
  dimnames(attr(mod, "diagnoses")) <- list(cd4stage = hDSlbl, agegr = hAGlbl, sex = sexlbl, year = yrlbl)
  dimnames(attr(mod, "late_diagnoses")) <- list(cd4stage = hDSlbl, agegr = hAGlbl, sex = sexlbl, year = yrlbl)
  dimnames(attr(mod, "artinit")) <- list(cd4stage = hDSlbl, agegr = hAGlbl, sex = sexlbl, year = yrlbl)
  dimnames(attr(mod, "artpop")) <- list(artdur = hTSlbl, cd4stage = hDSlbl,
                                        agegr = hAGlbl, sex = sexlbl, year = yrlbl)
  dimnames(attr(mod, "infections")) <- dimnames(mod)[-3]
  dimnames(attr(mod, "hivdeaths")) <- dimnames(mod)[-3]
  dimnames(attr(mod, "natdeaths")) <- dimnames(mod)[-3]
  names(attr(mod, "prev15to49")) <- yrlbl
  names(attr(mod, "incid15to49")) <- yrlbl
  mod
}

#' Simulate Shiny90 outputs saved in PJNZ file
#'
#' @param pjnz_path file path to .PJNZ file containing .shiny90
#'
#' @return Simulated model output from first90::simmod()
#' 
#' @details
#'
#' This function borrows code heavily from `naomi:::extract_shiny90_age_sex()`.
#'
#' It should extract output from a .shiny90 file saved with either R package
#' version or internal Spectrum Shiny90 implementation, which save the
#' components in slightly different formats.
#' 
shiny90_extract_sim <- function(pjnz_path) {

  print(basename(pjnz_path))

  shiny90_dir <- tempfile()
  on.exit(unlink(shiny90_dir))
  shiny90_name <- get_pjnz_shiny90_filename(pjnz_path)
  utils::unzip(pjnz_path, shiny90_name, exdir = shiny90_dir)
  shiny90_path <- file.path(shiny90_dir, shiny90_name)


  tmpd <- tempfile()
  on.exit(unlink(tmpd))
  utils::unzip(shiny90_path, exdir = tmpd)
  has_pjnz <- !is.null(pjnz_path) && file.exists(pjnz_path)

  if (file.exists(file.path(tmpd, "country.txt"))) {
    name <- brio::readLines(file.path(tmpd, "country.txt"))[1]
  } else if (has_pjnz) {
    name <- first90::read_country(pjnz_path)
    pjnz_region <- first90::read_region(pjnz_path)
    if (!is.null(pjnz_region)) {
      name <- paste0(name, " - ", pjnz_region)
    }
  } else {
    stop("PJNZ file required for .shiny90 created by Spectrum")
  }
  spectrum_data <- list.files(file.path(tmpd, "spectrum_data"), "rds$", full.names = TRUE)

  if (length(spectrum_data) > 0) {
    spec <- lapply(spectrum_data, readRDS)
    spec <- lapply(spec, "[[", "data")
  } else if (has_pjnz) {
    spec <- list(first90::extract_pjnz(pjnz_path))
  } else {
    stop("PJNZ file required for .shiny90 created by Spectrum")
  }

  fp <- first90::prepare_inputs_from_extracts(spec)
  if (!exists("popadjust", fp)) {
    fp$popadjust <- FALSE
  }

  if (file.exists(file.path(tmpd, "model_outputs/par.rds"))) {
    par <- readRDS(file.path(tmpd, "model_outputs/par.rds"))
  } else if (has_pjnz) {
    par <- as.numeric(readr::read_csv(file.path(tmpd, "model_outputs/par.csv"),
                                      col_names = FALSE, show_col_types = FALSE))
  } else {
    stop("PJNZ file required for .shiny90 created by Spectrum")
  }
  fpsim <- first90::create_hts_param(par, fp)
  mod <- first90::simmod(fpsim)

  mod <- first90_mod_add_dimnames(mod, fp$ss)
  attr(mod, "name") <- name

  mod
}

#' Format fine stratified Shiny90 simulation outputs in data frame
#'
#' @param mod model simulation output first90::simmod()
#'
#' @return
#' This returns a data frame with finely stratified Shiny90 simulation outputs.
#' Outputs are stratified by year, sex, coarse age group (15-16, 17-19, 20-24,
#' ..., 45-49, 50+), and CD4 stage.
#' 
#' Data frame contains outputs for:
#' 
#' * `hivpop_untreated`: Number of untreated PLHIV
#' * `artpop`: Number on ART (all durations; according to CD4 at initiation)
#' * `diagnpop_untreated`: Number of PLHIV who are diagnosed but untreated
#' * `diagnoses`: Number of diagnoses in the year
#' * `late_diagnoses`: Number of PLHIV who are drawn directly from undiagosed
#'      to on ART. These are counted in the number of `diagnoses` outcome.
#' * `artinit`: Number of ART initiations in the year
#'
shiny90_mod_extract_long <- function(mod) {
  modl <- attr(mod, "hivpop") %>%
    as.data.frame.table(responseName = "hivpop_untreated") %>%
    as_tibble() %>%
    left_join(
      attr(mod, "artpop") %>%
        as.data.frame.table(responseName = "artpop") %>%
        count(cd4stage, agegr, sex, year, wt = artpop, name = "artpop"),
      by = join_by(cd4stage, agegr, sex, year)
    ) %>%
    left_join(
      attr(mod, "diagnpop") %>%
        as.data.frame.table(responseName = "diagnpop_untreated"),
      by = join_by(cd4stage, agegr, sex, year)
    ) %>% 
    left_join(
      attr(mod, "diagnoses") %>%
        as.data.frame.table(responseName = "diagnoses"),
      by = join_by(cd4stage, agegr, sex, year)
    ) %>%   
    left_join(
      attr(mod, "late_diagnoses") %>%
        as.data.frame.table(responseName = "late_diagnoses"),
      by = join_by(cd4stage, agegr, sex, year)
    ) %>%
    left_join(
      attr(mod, "artinit") %>%
        as.data.frame.table(responseName = "artinit"),
      by = join_by(cd4stage, agegr, sex, year)
    )

  modl$name <- attr(mod, "name")

  modl <- modl %>%
    mutate(
      name = name
    ) %>%
    type.convert(as.is = TRUE) %>%
    select(name, year, sex, agegr, cd4stage, everything())

  modl
}
  

#' Aggregate outputs for CM model inputs
#'
#' @param modl data frame returned by
#'
#' @return
#' A data frame with aggregate for age 15+ with total
#' by {undiagnosed, untreated, on ART} and number with
#' CD4 <200 in each category.
#' 
aggregate_for_cm <- function(modl) {

  modl %>%
    mutate(
      cd4bel200 = as.integer(cd4stage %in% c("<50", "50-99", "100-199"))
    ) %>%
    summarise(
      undiagnosed = sum(hivpop_untreated) - sum(diagnpop_untreated),
      undiagnosed_cd4bel200 = sum(hivpop_untreated * cd4bel200) - sum(diagnpop_untreated * cd4bel200),
      diagnosed_untreated = sum(diagnpop_untreated),
      diagnosed_untreated_cd4bel200 = sum(diagnpop_untreated * cd4bel200),
      on_art = sum(artpop),
      diagnoses = sum(diagnoses),
      diagnoses_cd4bel200 = sum(diagnoses * cd4bel200),
      art_initiations = sum(artinit),
      art_initiations_cd4bel200 = sum(artinit * cd4bel200),
      .by = c(name, year)
    )

}


#' # Extract results
#'
#'

pjnz_dir <- "~/Data/Spectrum files/2023 final shared/SSA"
mwi_pjnz <- file.path(pjnz_dir, "Malawi_2023_National_HIV_estimates_Spectrum_AIM_model.pjnz")
uga_pjnz <- file.path(pjnz_dir, "Uganda-spectrum_20Feb2023_Calibratetosurvey.pjnz")
bwa_pjnz <- file.path(pjnz_dir, "Botswana2023v4 WPP 02_03_2023 KOS.PJNZ")
zwe_pjnz <- list.files(pjnz_dir, "^ZW_", full.names = TRUE)

out_for_cm <- c(mwi_pjnz, uga_pjnz, bwa_pjnz, zwe_pjnz) %>%
  lapply(shiny90_extract_sim) %>%
  lapply(shiny90_mod_extract_long) %>%
  lapply(aggregate_for_cm) %>%
  bind_rows() %>%
  filter(year %in% 2010:2022)

out_for_cm %>%
  filter(year == 2022) %>%
  mutate(plhiv_15plus = undiagnosed + diagnosed_untreated + on_art) %>%
  select(name, year, plhiv_15plus)

out_for_cm_nat <- out_for_cm %>%
  mutate(
    country = sub("(.*) - (.*)", "\\1", name)
  ) %>%
  summarise(across(-name, sum), .by = c(country, year))

out_for_cm_nat %>%
  filter(year == 2022) %>%
  mutate(plhiv_15plus = undiagnosed + diagnosed_untreated + on_art) %>%
  select(-contains("cd4bel200"))


out_for_cm_nat %>%
  filter(year %in% 2020:2022) %>%
  mutate(plhiv_15plus = undiagnosed + diagnosed_untreated + on_art) %>%
  select(-contains("cd4bel200"))

write_csv(out_for_cm_nat, "spectrum-2023_shiny90-untreated-below-200-stratification_2024-04-28.csv")

library(tidyr)
library(ggplot2)
library(scales)
library(forcats)
      
out_for_cm_nat %>%
  mutate(
    undiagnosed_propbel200 = undiagnosed_cd4bel200 / undiagnosed,
    diagnosed_untreated_propbel200 = diagnosed_untreated_cd4bel200 / diagnosed_untreated
  ) %>%
  pivot_longer(cols = c(undiagnosed_propbel200, diagnosed_untreated_propbel200)) %>%
  mutate(
    name = fct_relevel(name, "undiagnosed_propbel200", "diagnosed_untreated_propbel200")
  ) %>%
  ggplot(aes(year, value, color = country)) +
  geom_line() +
  scale_y_continuous("Percentage with CD4 <200", label = label_percent()) +
  scale_x_continuous(element_blank(), breaks = seq(2010, 2022, 2)) +
  facet_wrap(~name)

ggsave("undiagnosed-untreated-prop-below-cd4200_2024-04-28.pdf", w = 7.5, h = 3.5)