diff --git a/.gitignore b/.gitignore index cd724ac7..842374bb 100644 --- a/.gitignore +++ b/.gitignore @@ -21,3 +21,5 @@ tests/testthat/testdata/fit.RDS .vscode .Rprofile .idea + +.DS_Store diff --git a/DESCRIPTION b/DESCRIPTION index 455c4adc..c3341492 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: naomi Title: Naomi Model for Subnational HIV Estimates -Version: 2.9.23 +Version: 2.9.24 Authors@R: person(given = "Jeff", family = "Eaton", @@ -36,6 +36,7 @@ Imports: knitr, magrittr, mvtnorm, + naomi.resources (>= 0.0.3), naomi.options (>= 1.2.0), openxlsx, plotly, @@ -51,7 +52,6 @@ Imports: traduire, utils, withr, - writexl, yaml, zip, zoo @@ -78,8 +78,9 @@ Remotes: first90=mrc-ide/first90release, reside-ic/traduire, mrc-ide/naomi.options, + mrc-ide/naomi.resources, mrc-ide/mockr, mrc-ide/testthat.buildkite, - duckdb/duckdb-r@v0.9.1 + duckdb=duckdb/duckdb-r@v0.9.1 Config/testthat/edition: 3 Config/testthat/parallel: true diff --git a/NAMESPACE b/NAMESPACE index c7d328cb..0aba3158 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -44,9 +44,9 @@ export(get_plotting_metadata) export(hintr_calibrate) export(hintr_calibrate_plot) export(hintr_comparison_plot) -export(hintr_prepare_agyw_download) export(hintr_prepare_coarse_age_group_download) export(hintr_prepare_comparison_report_download) +export(hintr_prepare_shipp_download) export(hintr_prepare_spectrum_download) export(hintr_prepare_summary_report_download) export(hintr_run_model) diff --git a/NEWS.md b/NEWS.md index 8ac89764..4dd14b00 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,8 @@ +# naomi 2.9.24 + +* Generate PSE workbook from naomi outputs. +* Add function `extract_kp_workbook()` to extract key population totals from the Spectrum PJNZ file. Data are extracted from summary table saved in AIM Programme Statistics input for key populations. + # naomi 2.9.23 * Update Datim UIDs for Ethiopia 2024 boundary division. diff --git a/R/downloads.R b/R/downloads.R index 6216af85..5018cada 100644 --- a/R/downloads.R +++ b/R/downloads.R @@ -101,7 +101,7 @@ hintr_prepare_comparison_report_download <- function(output, ) } -#' Prepare AGYW tool download +#' Prepare SHIPP tool download #' #' @param hintr_output object #' @param path Path to save output file @@ -109,23 +109,31 @@ hintr_prepare_comparison_report_download <- function(output, #' #' @return Path to output file and metadata for file #' @export -hintr_prepare_agyw_download <- function(output, pjnz, +hintr_prepare_shipp_download <- function(output, pjnz, path = tempfile(fileext = ".xlsx")) { ## TODO: Do we need a version restriction on this? assert_model_output_version(output, "2.7.16") progress <- new_simple_progress() - progress$update_progress("PROGRESS_DOWNLOAD_AGYW") - dummy_data <- data.frame(x = c(1, 2, 3), y = c(3, 4, 5)) - writexl::write_xlsx(list(sheet = dummy_data), path = path) + progress$update_progress("PROGRESS_DOWNLOAD_SHIPP") + + risk_populations <- shipp_generate_risk_populations(output$model_output_path, + pjnz) + + sheets <- list( + "All outputs - F" = risk_populations$female_incidence, + "All outputs - M" = risk_populations$male_incidence, + "NAOMI outputs" = risk_populations$naomi_output + ) + write_shipp_workbook(sheets, dest = path) model_output <- read_hintr_output(output$model_output_path) options <- yaml::read_yaml(text = model_output$info$options.yml) list( path = path, metadata = list( - description = build_agyw_tool_description(options), + description = build_shipp_tool_description(options), areas = options$area_scope, - type = "agyw" + type = "shipp" ) ) } @@ -142,8 +150,8 @@ build_comparison_report_description <- function(options) { build_description(t_("DOWNLOAD_COMPARISON_DESCRIPTION"), options) } -build_agyw_tool_description <- function(options) { - build_description(t_("DOWNLOAD_AGYW_DESCRIPTION"), options) +build_shipp_tool_description <- function(options) { + build_description(t_("DOWNLOAD_SHIPP_DESCRIPTION"), options) } build_description <- function(type_text, options) { diff --git a/R/inputs-spectrum.R b/R/inputs-spectrum.R index c4107e1e..8af893e9 100644 --- a/R/inputs-spectrum.R +++ b/R/inputs-spectrum.R @@ -127,6 +127,8 @@ extract_pjnz_one <- function(pjnz, extract_shiny90) { #' @export #' extract_pjnz_program_data <- function(pjnz_list) { + + pjnz_list <- unroll_pjnz(pjnz_list) region_code <- lapply(pjnz_list, read_spectrum_region_code) @@ -185,7 +187,7 @@ read_dp_art_dec31 <- function(dp) { ## In Spectrum 2023, "" was updated to include children in the totals ## -> now need to sum over 5-year age groups for age 15+ to get the adult ART need - + male_15plus_needart <- dpsub("", 4:17*3 + 3, timedat.idx) male_15plus_needart <- vapply(lapply(male_15plus_needart, as.numeric), sum, numeric(1)) @@ -195,14 +197,14 @@ read_dp_art_dec31 <- function(dp) { art15plus_need <- rbind(male_15plus_needart, female_15plus_needart) dimnames(art15plus_need) <- list(sex = c("male", "female"), year = proj.years) - + if (any(art15plus_num[art15plus_isperc == 1] < 0 | art15plus_num[art15plus_isperc == 1] > 100)) { stop("Invalid percentage on ART entered for adult ART") } ## # Adult on ART adjustment factor - ## + ## ## * Implemented from around Spectrum 6.2 (a few versions before) ## * Allows user to specify scalar to reduce number on ART in each year ("") ## * Enabled / disabled by checkbox flag ("") @@ -281,7 +283,7 @@ read_dp_art_dec31 <- function(dp) { names(child_art) <- proj.years ## # Child on ART adjustment factor - ## + ## ## * Implemented same as adult adjustment factor above if (exists_dptag("") && @@ -684,7 +686,111 @@ extract_eppasm_pregprev <- function(mod, fp, years = NULL) { df } + + read_dp <- function(pjnz) { dpfile <- grep(".DP$", utils::unzip(pjnz, list = TRUE)$Name, value = TRUE) utils::read.csv(unz(pjnz, dpfile), as.is = TRUE) } + +#' Read key population summary data from PJNZ +#' +#' Reads key population summary data from Spectrum PJNZ. +#' +#' @param pjnz_list path to PJNZ file or zip of multiple PJNZ files +#' +#' +#' @noRd +#' + +extract_kp_workbook <- function(pjnz_list){ + + # Extract spectrum files + pjnz_list <- unroll_pjnz(pjnz_list) + + # Extract .DP files + dp <- lapply(pjnz_list, read_dp) + + # Extract kp workbook summary + kp <- lapply(dp, read_dp_keypop_summary) + + # Filter for spectrum file with consensus estimates + kp_out <- kp %>% + dplyr::bind_rows() %>% + dplyr::filter(!is.na(year)) + + # If no consensus estimates present, return empty dataframe + if(nrow(kp_out) == 0){kp_out <- kp[[1]]} + + # If multiple pjnz files, aggreagte consensus estimates + if(nrow(kp_out) > 4){ + + kp_out <- kp_out |> + dplyr::group_by(key_population, year, workbook_file) |> + dplyr::summarise(population_size = sum(population_size), + hiv_prevalence = sum(hiv_prevalence), + art_coverage = sum(art_coverage), + infections = sum(infections)) + + } + + kp_out + +} + + +#' Read key population summary data from PJNZ +#' +#' Reads key population summary data from Spectrum PJNZ. +#' +#' @param pjnz path to PJNZ file +#' +#' @examples +#' pjnz <- system.file("extdata/demo_mwi2019.PJNZ", package = "naomi") +#' dp <- dp <- naomi:::read_dp(pjnz) +#' read_dp_keypop_summary(dp) +#' +#' @noRd +#' +read_dp_keypop_summary <- function(dp) { + + exists_dptag <- function(tag, tagcol = 1) { + tag %in% dp[, tagcol] + } + dpsub <- function(tag, rows, cols, tagcol = 1) { + dp[which(dp[, tagcol] == tag) + rows, cols] + } + + kp_name <- c("FSW", "MSM", "TG", "PWID") + + if (exists_dptag("")) { + kp_tab <- dpsub("", 2:5, 4:7) + kp_tab <- sapply(kp_tab, as.numeric) + } else { + kp_tab <- matrix(NA, 4, 4) + } + + if (exists_dptag("")) { + kp_year <- as.integer(dpsub("", 2, 4)) + } else { + kp_year <- NA_integer_ + } + + if (exists_dptag("")) { + kp_file <- as.character(dpsub("", 2, 4)) + } else { + kp_file <- NA_character_ + } + + kp_summary <- data.frame( + key_population = kp_name, + year = kp_year, + population_size = kp_tab[1, ], + hiv_prevalence = kp_tab[2, ], + art_coverage = kp_tab[3, ], + infections = kp_tab[4, ], + workbook_file = kp_file + ) + + kp_summary +} diff --git a/R/shipp.R b/R/shipp.R new file mode 100644 index 00000000..c1494af0 --- /dev/null +++ b/R/shipp.R @@ -0,0 +1,1764 @@ +#' Format naomi outputs for PSE tool +#' +#' @param outputs Naomi output +#' @param options Naomi model options. +#' +#' +#' @return Naomi indicators formatted for the SHIPP workbook. +#' @keywords internal + + +shipp_format_naomi <- function(outputs, options){ + + naomi_ind <- outputs$indicators %>% + dplyr::filter(indicator %in% c("population", "plhiv", "infections","incidence", + "prevalence"), + calendar_quarter == options$calendar_quarter_t2) + + area_labels <- outputs$meta_area %>% + dplyr::select(area_id, area_name,area_level, spectrum_region_code) + area_label_cols <- intersect(names(naomi_ind), names(area_labels)) + + age_labels <- outputs$meta_age_group + age_label_cols <- intersect(names(naomi_ind), names(age_labels)) + + naomi_ind_labelled <- naomi_ind %>% + dplyr::left_join(area_labels, by = area_label_cols) %>% + dplyr::left_join(age_labels, by = age_label_cols) + + + summarise_naomi_ind <- function(dat, age_cat) { + + if(age_cat == "Y015_024"){age_groups <- c("Y015_019", "Y020_024")} + if(age_cat == "Y025_049"){age_groups <- c("Y025_029","Y030_034","Y035_039", + "Y040_044", "Y045_049")} + + dat %>% + dplyr::select(area_id, area_name, area_level, spectrum_region_code, calendar_quarter, + age_group, sex, indicator, mean) %>% + tidyr::pivot_wider(names_from = indicator, values_from = mean) %>% + dplyr::group_by(area_id, area_name, area_level,spectrum_region_code, + calendar_quarter, sex) %>% + dplyr::summarise( + "population" = sum(population * as.integer(age_group %in% age_groups)), + "plhiv" = sum(plhiv * as.integer(age_group %in% age_groups)), + "infections" = sum(infections * as.integer(age_group %in% age_groups)), + .groups = "drop") %>% + dplyr::mutate(age_group = age_cat, + incidence = (infections/(population - plhiv)), + prevalence = plhiv/population) %>% + tidyr::pivot_longer(cols = c(population, plhiv, infections, incidence, prevalence), + names_to = "indicator", + values_to = "mean") %>% + dplyr::mutate(age_group_label = dplyr::if_else(age_group == "Y015_024", "15-24", "25-49")) + } + + # Naomi indicators for aggregate age groups + df1 <- dplyr::bind_rows(summarise_naomi_ind(naomi_ind_labelled, "Y015_024"), + summarise_naomi_ind(naomi_ind_labelled, "Y025_049")) + + # Naomi indicators for 5-year age groups + 15-49 + df2 <- naomi_ind_labelled %>% + dplyr::filter(age_group %in% c("Y015_019", "Y020_024", "Y025_029", "Y030_034", + "Y035_039", "Y040_044", "Y045_049", "Y015_049")) %>% + dplyr::select(names(df1)) %>% + # Add aggregate indicators + dplyr::bind_rows(df1) %>% + # Format incidence from 1000 person years to 100 person years + dplyr::mutate(mean = dplyr::if_else(indicator == "incidence", mean * 100, mean)) + + # Format for workbook + df3 <- df2 %>% + dplyr::mutate(indicator = dplyr::recode(indicator, + "population" = "Pop", "plhiv" = "PLHIV", + "infections" = "new","incidence" = "Inci"), + sex = dplyr::recode(sex, + "female" = "f", "male" = "m", + "both" = "all"), + mean = as.character(mean)) + + # Incidence categories + df4 <- df3 %>% + dplyr::filter(indicator == "Inci") %>% + dplyr::mutate(mean = dplyr::case_when(mean < 0.3 ~ "Low", + mean >= 0.3 & mean< 1 ~ "Moderate", + mean >= 1 & mean < 3 ~ "High", + mean >= 3 ~ "Very High", + TRUE ~ NA_character_), + indicator = "Incicategory") + + # New infections for all age groups + sexes + df5 <- naomi_ind_labelled %>% + dplyr::filter(indicator == "infections", age_group == "Y000_999", sex == "both", + area_level == options$area_level) + + # Not all country names in meta_area match the country names in the spreadsheet - need + # to match to populate tabs + + country_name_db <- tibble::tribble(~country, ~iso3, + "Botswana", "BWA", + "Cameroon", "CMR", + "Kenya", "KEN", + "Lesotho", "LSO", + "Mozambique", "MOZ", + "Malawi", "MWI_demo", + "Malawi", "MWI", + "Namibia", "NAM", + "Eswatini", "SWZ", + "Tanzania", "TZA", + "Uganda", "UGA", + "South Africa", "ZAF", + "Zambia", "ZMB", + "Zimbabwe", "ZWE", + "Angola", "AGO", + "Burundi", "BDI", + "Democratic Republic of the Congo", "COD", + "Gabon", "GAB", + "Rwanda", "RWA", + "Ethiopia", "ETH", + "Haiti", "HTI", + "Chad", "TCD", + "Cote D'Ivoire", "CIV", + "Ghana", "GHA", + "Guinea", "GIN", + "Liberia", "LBR", + "Mali", "MLI", + "Sierra Leone", "SLE", + "Togo", "TGO", + "Burkina Faso", "BFA", + "Congo", "COG", + "Benin", "BEN", + "Central African Republic", "CAF", + "The Gambia", "GMB", + "Guinea bissau", "GNB", + "Equatorial Guinea", "GNQ", + "Niger", "NER", + "Nigeria", "NGA", + "Senegal", "SEN") + + country_name <- country_name_db[country_name_db$iso3 == options$area_scope,]$country + + # Format + naomi_wide <- dplyr::bind_rows(df3, df4) %>% + dplyr::filter(area_level == options$area_level) %>% + tidyr::pivot_wider(id_cols = c(area_id,area_name), + names_from = c(indicator,age_group_label,sex), + names_sep = "", values_from = mean) %>% + dplyr::mutate(Country = country_name, newAll = df5$mean) %>% + dplyr::select(Country,area_id,area_name,`Pop15-24all`,`Pop15-24f`,`Pop15-24m`, + `PLHIV15-24all`,`PLHIV15-24f`,`PLHIV15-24m`, + newAll, `new15-24all`,`new15-24f`,`new15-24m`, + `Inci15-24f`,`Incicategory15-24f`,`Inci15-24m`,`Incicategory15-24m`, + `Pop15-19all`,`Pop15-19f`,`Pop15-19m`, + `PLHIV15-19all`,`PLHIV15-19f`,`PLHIV15-19m`, + `new15-19all`,`new15-19f`,`new15-19m`, + `Inci15-19f`,`Incicategory15-19f`,`Inci15-19m`,`Incicategory15-19m`, + `Pop20-24all`,`Pop20-24f`,`Pop20-24m`, + `PLHIV20-24all`,`PLHIV20-24f`,`PLHIV20-24m`, + `new20-24all`,`new20-24f`,`new20-24m`, + `Inci20-24f`,`Incicategory20-24f`,`Inci20-24m`,`Incicategory20-24m`, + `Pop25-49all`,`Pop25-49f`,`Pop25-49m`, + `PLHIV25-49all`,`PLHIV25-49f`,`PLHIV25-49m`, + `new25-49all`,`new25-49f`,`new25-49m`, + `Inci25-49f`,`Incicategory25-49f`,`Inci25-49m`,`Incicategory25-49m`, + `Pop25-29all`,`Pop25-29f`,`Pop25-29m`, + `PLHIV25-29all`,`PLHIV25-29f`,`PLHIV25-29m`, + `new25-29all`,`new25-29f`,`new25-29m`, + `Inci25-29f`,`Incicategory25-29f`,`Inci25-29m`,`Incicategory25-29m`, + `Pop30-34all`,`Pop30-34f`,`Pop30-34m`, + `PLHIV30-34all`,`PLHIV30-34f`,`PLHIV30-34m`, + `new30-34all`,`new30-34f`,`new30-34m`, + `Inci30-34f`,`Incicategory30-34f`,`Inci30-34m`,`Incicategory30-34m`, + `Pop35-39all`,`Pop35-39f`,`Pop35-39m`, + `PLHIV35-39all`,`PLHIV35-39f`,`PLHIV35-39m`, + `new35-39all`,`new35-39f`,`new35-39m`, + `Inci35-39f`,`Incicategory35-39f`,`Inci35-39m`,`Incicategory35-39m`, + `Pop40-44all`,`Pop40-44f`,`Pop40-44m`, + `PLHIV40-44all`,`PLHIV40-44f`,`PLHIV40-44m`, + `new40-44all`,`new40-44f`,`new40-44m`, + `Inci40-44f`,`Incicategory40-44f`,`Inci40-44m`,`Incicategory40-44m`, + `Pop45-49all`,`Pop45-49f`,`Pop45-49m`, + `PLHIV45-49all`,`PLHIV45-49f`,`PLHIV45-49m`, + `new45-49all`,`new45-49f`,`new45-49m`, + `Inci45-49f`,`Incicategory45-49f`,`Inci45-49m`,`Incicategory45-49m`, + `Pop15-49all`,`Pop15-49f`,`Pop15-49m`, + `PLHIV15-49all`,`PLHIV15-49f`,`PLHIV15-49m`, + `new15-49all`,`new15-49f`,`new15-49m`, + `Inci15-49f`,`Incicategory15-49f`,`Inci15-49m`,`Incicategory15-49m`) + + ## Clean up area names and Country names + if(options$area_scope=="AGO") { + naomi_wide$area_name <- stringr::str_to_title(naomi_wide$area_name) + } + + if(options$area_scope %in% c("TCD","GIN")) { + naomi_wide$area_name <- iconv(naomi_wide$area_name, from="UTF-8",to="LATIN1") + } + + v <- list(naomi_long = df2, + naomi_wide = naomi_wide) + + } + + + +#' Dissagreggate admin1 FSW proportions from Oli's KP model to 5-age groups +#' +#' @param outputs Naomi output. +#' @param options Naomi model options. +#' @param naomi_pop Naomi population estimates for T2. +#' @param kp_consensus Key pop consensus estimates. +#' +#' +#' @return District level FSW estimates by 5-year age bands for ages 15-49. +#' @keywords internal + + +shipp_disaggregate_fsw <- function(outputs, + options, + naomi_pop, + kp_consensus){ + + # Extract country specific national FSW PSEs + iso3 <- options$area_scope + + pse <- naomi.resources::load_shipp_exdata("kp_estimates", iso3) %>% + dplyr::filter(kp == "FSW", indicator == "pse_prop") + + fsw_pse <- pse %>% + dplyr::rename(prop_fsw = median) %>% + dplyr::select(-indicator,-lower,-upper) + + age_groups <- c("Y015_019", "Y020_024", "Y025_029", "Y030_034", + "Y035_039", "Y040_044", "Y045_049") + + # Calculating FSW proportion of total female population + fsw <- fsw_pse %>% + dplyr::mutate(age_group = "Y015_049") %>% + dplyr::left_join(naomi_pop %>% dplyr::filter(sex == "female"), + by = dplyr::join_by(iso3, area_id, age_group)) %>% + dplyr::mutate(total_fsw = population * prop_fsw) %>% + dplyr::select(iso3, area_id, total_fsw, age_group, area_level, spectrum_region_code) + + # Check for consensus estimate of FSW + fsw_consensus <- kp_consensus[kp_consensus$key_population == "FSW", ]$population_size + + if(!is.na(fsw_consensus) & fsw_consensus > 0){ + + # Check if consensus estimate is larger than age matched population denominator + pop <- naomi_pop[naomi_pop$area_level == 0 & naomi_pop$age_group == "Y015_049" & naomi_pop$sex == "female",]$population + prop_pop <- fsw_consensus / pop + + if(prop_pop < 0.05) { + + # Scale total FSW population to consensus PSE estimate + fsw_scaled <- fsw %>% + dplyr::mutate( + relative_prop = total_fsw/sum(total_fsw), + consensus_pse = fsw_consensus, + total_fsw = consensus_pse * relative_prop) + + fsw <- fsw_scaled %>% dplyr::select(-consensus_pse, relative_prop) + } + + } + + + # FSW age distribution parameters in ZAF from Thembisa + # Downloaded from: https://www.thembisa.org/content/downloadPage/Thembisa4_3 + gamma_mean <- 29 + gamma_sd <- 9 + beta <- gamma_mean / gamma_sd^2 # rate + alpha <- gamma_mean * beta # shape + + # Distribution function of the gamma + zaf_gamma <- data.frame( + dist = diff(pgamma(c(15, 20, 25, 30, 35, 40, 45, 50), shape = alpha, rate = beta)), + age_group = age_groups) %>% + dplyr::mutate(dist = dist / sum(dist)) + + pskewlogis <- function(t, scale, shape, skew) { + (1 + (scale * t)^-shape)^-skew + } + + # Calculate proportion of sexually active population using Kinh's country specific + # estimates of age at first sex and naomi population + afs <- naomi.resources::load_shipp_exdata("afs", iso3) + + # Select birth cohort from 2000, to turn 15 in 2015 + cohort <- 2000 + + afs <- afs %>% + dplyr::filter(yob == cohort, sex == "female", iso3 == options$area_scope) %>% + dplyr::full_join(dplyr::select(fsw,iso3,area_id), multiple = "all", by = dplyr::join_by(iso3)) + + df <- data.frame() + + # Calculate sexually active population by age and sex for each district + for (x in unique(afs$area_id)) { + afs_x <- dplyr::filter(afs, area_id == x) + ages <- 15:49 + + df_x <- data.frame( + area_id = x, + age = ages, + eversex = pskewlogis( + ages, + scale = afs_x$lambda, + skew = afs_x$skew, + shape = afs_x$shape + ), + age_group = rep(age_groups, each = 5) + ) + + df_x <- df_x %>% + dplyr::group_by(area_id, age_group) %>% + dplyr::summarise(eversex = mean(eversex), .groups = "drop") %>% + dplyr::left_join( + naomi_pop %>% dplyr::filter(sex == "female"), + by = c("area_id", "age_group") + ) %>% + dplyr::mutate( + eversexpop = eversex * population, + eversexpop_prop = eversexpop / sum(eversexpop) + ) + + df <- dplyr::bind_rows(df, df_x) + } + + # Adjusting country specific sexual debut estimates with age distribution of + # FSW from Thembisa + #Downloaded from: https://www.thembisa.org/content/downloadPage/Thembisa4_3 + zaf_propensity <- naomi.resources::load_shipp_exdata("zaf_propensity", iso3 = "ZAF") %>% + dplyr::filter(kp == "FSW") + + fsw_est <- df %>% + # Add FSW propensity estimates from ZAF + dplyr::left_join(zaf_propensity, by = "age_group") %>% + # Calculate distribution of FSWs + dplyr::mutate(dist = eversexpop_prop * propensity) %>% + dplyr::group_by(area_id) %>% + dplyr::mutate(dist = dist / sum(dist)) %>% + dplyr::ungroup() %>% + # Add FSW PSEs + dplyr::full_join( + fsw %>% dplyr::select(total_fsw, iso3, area_id, area_level), + by = dplyr::join_by(area_id, iso3, area_level) + ) %>% + # Calculate FSW proportions + dplyr::mutate( + fsw = dist * total_fsw, + fsw_prop = fsw / population, + consensus_estimate = fsw_consensus + ) %>% + dplyr::select(-eversexpop, -eversexpop_prop, -propensity, -dist, -total_fsw) + + fsw_est + +} + + +#' Disaggregate admin1 PWID proportions from Oli's KP model to 5-age groups +#' +#' @param outputs Naomi output. +#' @param options Naomi model options. +#' @param naomi_pop Naomi population estimates for T2. +#' @param kp_consensus Key pop consensus estimates. +#' +#' @return District level PWID estimates by 5-year age bands for ages 15-49. +#' @keywords internal + +shipp_disaggregate_pwid <- function(outputs, + options, + naomi_pop, + kp_consensus) +{ + + # Extract country specific national PWID PSEs + iso3 <- options$area_scope + + pse <- naomi.resources::load_shipp_exdata("kp_estimates", iso3) %>% + dplyr::filter(kp == "PWID", indicator == "pse_prop") + + pwid_pse <- pse %>% + dplyr::rename(prop_pwid = median) %>% + dplyr::select(-indicator,-lower,-upper) + + age_groups <- c("Y015_019", "Y020_024", "Y025_029", "Y030_034", + "Y035_039", "Y040_044", "Y045_049") + + pwid <- pwid_pse %>% + dplyr::mutate(age_group = "Y015_049") %>% + dplyr::left_join(naomi_pop %>% dplyr::filter(sex == "male"), + by = dplyr::join_by(iso3, area_id, age_group)) %>% + dplyr::mutate(total_pwid = population * prop_pwid) %>% + dplyr::select(iso3, area_id, total_pwid, age_group, area_level) + + # Check for consensus estimate of PWID + pwid_consensus <- kp_consensus[kp_consensus$key_population == "PWID", ]$population_size + + if(!is.na(pwid_consensus) & pwid_consensus > 0){ + # Check if consensus estimate is larger than age matched population denominator + pop <- naomi_pop[naomi_pop$area_level == 0 & naomi_pop$age_group == "Y015_049" & naomi_pop$sex == "male",]$population + prop_pop <- pwid_consensus / pop + + if(prop_pop < 0.05) { + + # Scale total PWID population to consensus PSE estimate + pwid_scaled <- pwid %>% + dplyr::mutate( + relative_prop = total_pwid/sum(total_pwid), + consensus_pse = pwid_consensus, + total_pwid = consensus_pse * relative_prop) + + pwid <- pwid_scaled %>% dplyr::select(-consensus_pse, relative_prop) + } + + } + + + + # Assumption from literature that 9% of PWID are female so remove them from + # the male denominator + + pwid$total_pwid <- pwid$total_pwid * 0.91 + + # PWID age distribution + # Review of literature - Hines et al Lancet Global Health 2020: + # DOI:https://doi.org/10.1016/S2214-109X(19)30462-0 + gamma_mean <- 29.4 + gamma_sd <- 7 + beta <- gamma_mean / gamma_sd^2 # rate + alpha <- gamma_mean * beta # shape + + # Distribution function of the gamma + zaf_gamma <- data.frame( + dist = diff(pgamma(c(15, 20, 25, 30, 35, 40, 45, 50), shape = alpha, rate = beta)), + age_group = age_groups) %>% + dplyr::mutate(dist = dist / sum(dist)) + + + # Naomi population + pop <- naomi_pop %>% + dplyr::filter(area_id %in% unique(pwid$area_id), + age_group %in% age_groups, + sex == "male") + + pwid_est <- dplyr::left_join( + pop, zaf_gamma, + by = dplyr::join_by(age_group)) %>% + dplyr::full_join( + dplyr::select(pwid,total_pwid, iso3, area_id), + by = c("area_id", "iso3") + ) %>% + dplyr::mutate(pwid = dist * total_pwid, + pwid_prop = pwid / population, + consensus_estimate = pwid_consensus) %>% + dplyr::select( -dist, -total_pwid, -sex) + + pwid_est +} + +#' Disaggregate admin1 MSM proportions from Oli's KP model to 5-age groups +#' +#' @param outputs Naomi output. +#' @param options Naomi model options. +#' @param naomi_pop Naomi population estimates for T2. +#' @param kp_consensus Key pop consensus estimates. +#' +#' @return District level MSM estimates by 5-year age bands for ages 15-49. +#' @keywords internal + +shipp_disaggregate_msm <- function(outputs, + options, + naomi_pop, + kp_consensus) +{ + + # Extract country specific national MSM PSEs + iso3 <- options$area_scope + pse <- naomi.resources::load_shipp_exdata("kp_estimates", iso3) %>% + dplyr::filter(kp == "MSM", indicator == "pse_prop") + + msm_pse <- pse %>% + dplyr::rename(prop_msm = median) %>% + dplyr::select(-indicator,-lower,-upper) + + age_groups <- c("Y015_019", "Y020_024", "Y025_029", "Y030_034", + "Y035_039", "Y040_044", "Y045_049") + + msm <- msm_pse %>% + dplyr::mutate(age_group = "Y015_049") %>% + dplyr::left_join(naomi_pop %>% dplyr::filter(sex == "male"), + by = dplyr::join_by(iso3, area_id, age_group)) %>% + dplyr::mutate(total_msm = population * prop_msm) %>% + dplyr::select(iso3, area_id, total_msm, age_group, area_level) + + # Check for consensus estimate of MSM + msm_consensus <- kp_consensus[kp_consensus$key_population == "MSM", ]$population_size + + if(!is.na(msm_consensus) & msm_consensus > 0){ + + # Check if consensus estimate is larger than age matched population denominator + pop <- naomi_pop[naomi_pop$area_level == 0 & naomi_pop$age_group == "Y015_049" & naomi_pop$sex == "male",]$population + prop_pop <- msm_consensus / pop + + if(prop_pop < 0.05) { + + # Scale total MSM population to consensus PSE estimate + msm_scaled <- msm %>% + dplyr::mutate( + relative_prop = total_msm/sum(total_msm), + consensus_pse = msm_consensus, + total_msm = consensus_pse * relative_prop) + + msm <- msm_scaled %>% dplyr::select(-consensus_pse, relative_prop) + } + + } + + + + + + # MSM age distribution parameters in ZAF from Thembisa + # Downloaded from: https://www.thembisa.org/content/downloadPage/Thembisa4_3report + gamma_mean <- 25 + gamma_sd <- 7 + beta <- gamma_mean / gamma_sd^2 # rate + alpha <- gamma_mean * beta # shape + + # Distribution function of the gamma + zaf_gamma <- data.frame( + dist = diff(pgamma(c(15, 20, 25, 30, 35, 40, 45, 50), shape = alpha, rate = beta)), + age_group = age_groups) %>% + dplyr::mutate(dist = dist / sum(dist)) + + + pskewlogis <- function(t, scale, shape, skew) { + (1 + (scale * t)^-shape)^-skew + } + + # Calculate proportion of sexually active population using Kinh's country specific + # estimates of age at first sex and naomi population + afs <- naomi.resources::load_shipp_exdata("afs", iso3) + + # Select birth cohort from 2000, to turn 15 in 2015 + cohort <- 2000 + + afs <- afs %>% + dplyr::filter(yob == cohort, sex == "male", iso3 == options$area_scope) %>% + dplyr::full_join(dplyr::select(msm,iso3,area_id), multiple = "all", by = dplyr::join_by(iso3)) + + df <- data.frame() + + # Calculate sexually active population by age and sex for each district + for(x in unique(afs$area_id)) { + afs_x <- dplyr::filter(afs, area_id == x) + ages <- 15:49 + + df_x <- data.frame( + area_id = x, + age = ages, + eversex = pskewlogis( + ages, + scale = afs_x$lambda, + skew = afs_x$skew, + shape = afs_x$shape + ), + age_group = rep(age_groups, each = 5) + ) + + df_x <- df_x %>% + dplyr::group_by(area_id, age_group) %>% + dplyr::summarise(eversex = mean(eversex), .groups = "drop") %>% + dplyr::left_join( + naomi_pop %>% dplyr::filter(sex == "male"), + by = c("area_id", "age_group") + ) %>% + dplyr::mutate( + eversexpop = eversex * population, + eversexpop_prop = eversexpop / sum(eversexpop) + ) + + df <- dplyr::bind_rows(df, df_x) + } + + # Adjusting country specific sexual debut estimates with age distribution of + # MSM from Thembisa + zaf_propensity <- naomi.resources::load_shipp_exdata("zaf_propensity", iso3 = "ZAF") %>% + dplyr::filter(kp == "MSM") + + + msm_est <- df %>% + # Add MSM propensity estimates from ZAF + dplyr::left_join(zaf_propensity, by = "age_group") %>% + # Calculate distribution of MSM + dplyr::mutate(dist = eversexpop_prop * propensity) %>% + dplyr::group_by(area_id) %>% + dplyr::mutate(dist = dist / sum(dist)) %>% + dplyr::ungroup() %>% + # Add MSM PSEs + dplyr::full_join( + msm %>% dplyr::select(total_msm, iso3, area_id, area_level), + by = dplyr::join_by(area_id, iso3, area_level) + ) %>% + # Calculate FSW proportions + dplyr::mutate(msm = dist * total_msm, + msm_prop = msm / population, + consensus_estimate = msm_consensus) %>% + dplyr::select(-eversexpop, -eversexpop_prop, -propensity, - dist, -total_msm) + + msm_est +} + +#' Adjust female sexual behavior risk groups by FSW proportions +#' +#' @param outputs Naomi output. +#' @param options Naomi model options. +#' @param fsw_est 5-year estimates of FSW PSEs generated from `shipp_disaggregate_fsw()`. +#' +#' +#' Estimates are generated for the following groups: +#' +#' * `nosex12m`: +#' * `sexcohab`: +#' * `sexnonregplus`: +#' * `sexnonreg`: +#' * `sexpaid12m`: +#' * `nosex12m`: +#' +#' Calculation steps: +#' 1. Align admin0/admin1 FSW proportions with SRB SAE estimates. +#' 2. Subtract the proportion of FSW from total high risk female population. +#' +#' @return District level estimates of female sexual risk behaviour groups +#' @keywords internal + + +shipp_adjust_sexbehav_fsw <- function(outputs, + options, + fsw_est) +{ + + + # Match FSW estimates (admin0 or admin1) with SAE estimates + fsw_analysis_level <- paste0("area_id",unique(fsw_est$area_level)) + + areas_wide <- naomi::spread_areas(outputs$meta_area) %>% + sf::st_drop_geometry() + + map <- dplyr::select(areas_wide, area_id, dplyr::all_of(fsw_analysis_level)) %>% + dplyr::rename(fsw_match_area = 2) + + # Allocate admin1 FSW proportions + fsw_df <- fsw_est %>% dplyr::select(age_group, fsw_match_area = area_id, fsw_prop) + + # Load female SRB proportions + female_srb <- naomi.resources::load_shipp_exdata("srb_female", options$area_scope) + + adj_female_srb <- female_srb %>% + dplyr::filter(iso3 == options$area_scope) %>% + dplyr::left_join(map, by = dplyr::join_by(area_id)) %>% + dplyr::left_join(fsw_df, by = dplyr::join_by(age_group, fsw_match_area)) %>% + dplyr::select(iso3, year, indicator, survey_id, age_group, area_id, area_name, estimate_smoothed, fsw_prop) %>% + tidyr::pivot_wider( + names_from = indicator, + values_from = estimate_smoothed + ) %>% + # Subtracting proportion FSW from total high risk female population + dplyr::mutate( + sexpaid12m = fsw_prop, + sexnonreg = 1 - nosex12m - sexcohab - sexpaid12m, + fsw_prop = NULL + ) %>% + tidyr::pivot_longer( + cols = c(nosex12m,sexcohab,sexnonregplus,sexnonreg,sexpaid12m), + names_to = "indicator", + values_to = "estimate_smoothed" + ) + + adj_female_srb + +} + +#' Adjust male sexual behavior risk groups by MSM + PWID proportions +#' +#' @param outputs Naomi output. +#' @param options Naomi model options. +#' @param msm_est 5-year estimates of MSM PSEs generated from `shipp__disaggregate_msm()`. +#' @param pwid_est 5-year estimates of MSM PSEs generated from `shipp__disaggregate_pwid()`. +#' +#' @return District level estimates of male sexual risk behaviour groups +#' @keywords internal +#' +#' Estimates are generated for the following groups: +#' +#' * `nosex12m`: +#' * `sexcohab`: +#' * `sexonregplus`: +#' * `sexonreg`: +#' * `msm`: +#' * `pwid`: +#' +#' +#' Calculation steps: +#' 1. Align admin0/admin1 MSM and PWID proportions with SRB SAE estimates. +#' 2. Subtracting MSM and PWID proportionally from all SRB groups. + + +shipp_adjust_sexbehav_msm_pwid <- function(outputs, + options, + msm_est, + pwid_est) { + + + # Match KP estimates (admin0 or admin1) with SAE estimates + msm_analysis_level <- paste0("area_id",unique(msm_est$area_level)) + + areas_wide <- naomi::spread_areas(outputs$meta_area) %>% + sf::st_drop_geometry() + + map <- dplyr::select(areas_wide, area_id, dplyr::all_of(msm_analysis_level)) %>% + dplyr::rename(kp_match_area = 2) + + # Allocate KP + msm_df <- msm_est %>% dplyr::select(age_group, kp_match_area = area_id, msm_prop) + pwid_df <- pwid_est %>% dplyr::select(age_group, kp_match_area = area_id, pwid_prop) + + # Load male SRB proportions + male_srb <- naomi.resources::load_shipp_exdata("srb_male", options$area_scope) + + adj_male_srb <- male_srb %>% + dplyr::filter(iso3 == options$area_scope) %>% + dplyr::left_join(map, by = dplyr::join_by(area_id)) %>% + dplyr::left_join(msm_df, by = dplyr::join_by(age_group, kp_match_area)) %>% + dplyr::left_join(pwid_df, by = dplyr::join_by(age_group, kp_match_area)) %>% + dplyr::select(iso3, year, indicator, survey_id, age_group, area_id, area_name, + estimate_smoothed, msm_prop, pwid_prop) %>% + tidyr::pivot_wider( + names_from = indicator, + values_from = estimate_smoothed + ) %>% + # Subtracting MSM and PWID proportionally from all SRB risk groups + # (FSW was just from high-risk females) + dplyr::mutate( + nosex12m = nosex12m * (1 - pwid_prop - msm_prop), + sexcohab = sexcohab * (1 - pwid_prop - msm_prop), + sexnonregplus = sexnonregplus * (1 - pwid_prop - msm_prop), + sexnonreg = sexnonregplus, + msm = msm_prop, msm_prop = NULL, + pwid = pwid_prop, pwid_prop = NULL + ) %>% + tidyr::pivot_longer( + cols = c(nosex12m, sexcohab, sexnonregplus, sexnonreg, msm, pwid), + names_to = "indicator", + values_to = "estimate_smoothed" + ) + + adj_male_srb + +} + +#' Calculate prevalence for female SRB groups. +#' +#' @param naomi_output Naomi output. +#' @param options Naomi model options. +#' @param fsw_est 5-year estimates of FSW PSEs generated from `shipp_disaggregate_fse()`. +#' @param female_srb FSW adjusted estimates of female SRB groups generated by `shipp_adjust_sexbehav_fsw()` +#' @param survey_year Year of survey to sample estimates. +#' +#' To calculate district-age-sex-sexual behaviour-specific HIV prevalence, we maintain +#' HIV prevalence from Naomi for a district-age-sex, but disaggregate to different +#' risk behaviours using: +#' (1) HIV prevalence ratios from household surveys for those reporting no +#' sex vs one cohabiting vs non-regular sexual partner(s), and +#' (2) a linear regression through admin-1 level estimates of the ratio of KP +#' prevalence to gen-pop prevalence used to predict an age-district-specific +#' FSW to general population prevalence ratio. +#' +#' @return SRB PSEs with logit prevalence estimates. +#' @keywords internal + +shipp_calculate_prevalence_female <- function(naomi_output, + options, + fsw_est, + female_srb, + survey_year_sample = 2018) { + + # Naomi estimates of PLHIV and population by district and age band + naomi_est <- naomi_output %>% + dplyr::filter(calendar_quarter == options$calendar_quarter_t2, + sex == "female", + indicator %in% c("population", "plhiv", "infections", "prevalence")) %>% + dplyr::select(area_id, area_level, age_group, indicator, mean) %>% + tidyr::pivot_wider(names_from = indicator, values_from = mean) %>% + dplyr::rename(gen_prev = prevalence) + + # Naomi general population prevalence + genpop_prev <- naomi_est %>% + dplyr::filter(age_group == "Y015_049") %>% + dplyr::select(area_id, gen_prev) + + # Extract country specific national FSW prevalence + iso3 <- options$area_scope + # THIS IS NOW USING SINGLE COUNTRY INSTEAD OF ALL COUNTRIES + fsw_prev <- naomi.resources::load_shipp_exdata("kp_estimates", iso3) %>% + dplyr::filter(kp == "FSW", indicator == "prevalence") + + kp_prev <- fsw_prev %>% + dplyr::select(iso3,area_id,median) %>% + dplyr::left_join(genpop_prev, by = dplyr::join_by(area_id)) %>% + dplyr::mutate(prev_fsw_logodds = log(median / (1 - median)), + prev_logodds = log(gen_prev / (1 - gen_prev))) + + # KP regression: FSW prevalence relative to general prevalence + # ########## THIS REGRESSION SHOULD BE TAKING DATA FROM ALL ADMIN-1 LEVEL + kp_fit <- lm(prev_fsw_logodds ~ prev_logodds, data = kp_prev) + + # Modelled estimates of proportion in each risk group + risk_group_prop <- female_srb %>% + dplyr::filter(year == survey_year_sample) %>% + dplyr::select(area_id, age_group, indicator, estimate_smoothed) %>% + tidyr::pivot_wider( names_from = indicator, values_from = estimate_smoothed, + values_fn = mean, names_prefix = "prop_") %>% + dplyr::right_join(naomi_est, by = c("area_id", "age_group")) %>% + dplyr::filter(!is.na(prop_nosex12m)) %>% + dplyr::mutate( + population_nosex12m = population * prop_nosex12m, + population_sexcohab = population * prop_sexcohab, + population_sexnonreg = population * prop_sexnonreg, + population_sexpaid12m = population * prop_sexpaid12m + ) + + # Calculate prevalence in each category + calculate_prevalence <- function(x, iso3){ + + # Log odds ratio from SRB group survey prevalence + lor <- naomi.resources:::load_shipp_exdata("srb_survey_lor", iso3) %>% + dplyr::filter(sex == "female") + + lor_15to29 <- lor$lor_15to29 + names(lor_15to29) <- lor$srb_group + + lor_30to49 <- lor$lor_30to49 + names(lor_30to49) <- lor$srb_group + + if (x$age_group[1] %in% c("Y015_019","Y020_024","Y025_029")) { + lor <- lor_15to29 + } else { + lor <- lor_30to49 + } + + population_fine <- dplyr::filter(x, indicator == "population")$estimate + plhiv <- x$plhiv[1] + ywkp_lor <- c("ywkp_lor" = x$ywkp_lor[1]) + lor <- c(lor[1:3],ywkp_lor) + prev <- logit_scale_prev(lor, population_fine, plhiv) + y <- dplyr::filter(x, indicator == "prop") %>% + dplyr::mutate( indicator = "prev", estimate = prev) + dplyr::bind_rows(x, y) + } + + # Calculate logit prevalence and format + logit_prev <- risk_group_prop %>% + dplyr::mutate(ywkp_lor = calculate_ywkp_pr_lor(gen_prev, fit = kp_fit)$lor) %>% + dplyr::select(-starts_with("pr_"), -gen_prev) %>% + tidyr::pivot_longer( + cols = starts_with(c("population_", "prop_")), + names_to = "indicator", + values_to = "estimate" + ) %>% + tidyr::separate( + indicator, + into = c("indicator", "behav") + ) %>% + dplyr::filter(behav %in% c("nosex12m", "sexcohab", "sexnonreg", "sexpaid12m")) %>% + split(~ area_id + age_group) %>% + lapply(calculate_prevalence, iso3) %>% + dplyr::bind_rows() %>% + tidyr::unite("indicator", indicator, behav, sep = "_") %>% + tidyr::pivot_wider( names_from = indicator, values_from = estimate) %>% + dplyr::mutate_if(is.numeric, as.numeric) %>% + dplyr::mutate_if(is.factor, as.character) + + logit_prev + + +} + +#' Calculate prevalence for male SRB groups. +#' +#' @param naomi_output Naomi output. +#' @param areas Naomi boundary file. +#' @param options Naomi model options. +#' @param msm_est 5-year estimates of MSM PSEs generated from `shipp__disaggregate_msm()`. +#' @param male_srb MSM and PWID adjusted estimates of male SRB groups generated by `shipp_adjust_sexbehav_msm_pwid()`. +#' @param survey_year Year of survey to sample estimates. +#' +#' +#' To calculate district-age-sex-sexual behaviour-specific HIV prevalence, we maintain +#' HIV prevalence from Naomi for a district-age-sex, but disaggregate to different +#' risk behaviours using: +#' (1) HIV prevalence ratios from household surveys for those reporting no +#' sex vs one cohabiting vs non-regular sexual partner(s), and +#' (2) admin-1 level estimates of the ratio of KP prevalence to gen-pop prevalence +#' among 15-24 year olds for MSM (due to the young age distribution of MSM) or +#' among 15-49 year olds for PWID (due to the older age distribution of PWID) +#' applied to all age groups among MSM and PWID in districts by admin-1 unit. +#' +#' @return SRB PSEs with logit prevalence estimates. +#' @keywords internal + + +shipp_calculate_prevalence_male <- function(naomi_output, + areas, + options, + msm_est, + male_srb, + survey_year_sample = 2018) { + + # Naomi estimates of PLHIV and population by district and age band + naomi_est <- naomi_output %>% + dplyr::filter(calendar_quarter == options$calendar_quarter_t2, + sex == "male", + indicator %in% c("population", "plhiv", "infections", "prevalence")) %>% + dplyr::select(area_id, area_level, age_group, indicator, mean) %>% + tidyr::pivot_wider(names_from = indicator, values_from = mean) %>% + dplyr::rename(gen_prev = prevalence) + + # Naomi general population prevalence + genpop_prev <- naomi_est %>% + dplyr::filter(age_group == "Y015_024" | age_group == "Y015_049") %>% + dplyr::select(area_id, area_level, age_group, gen_prev) %>% + tidyr::pivot_wider(names_from = age_group, values_from = gen_prev) %>% + dplyr::mutate(logit_gen_prev_msm = log(Y015_024 / (1 - Y015_024)), + logit_gen_prev_pwid = log(Y015_049 / (1 - Y015_049))) %>% + dplyr::select(area_id, logit_gen_prev_msm, logit_gen_prev_pwid, area_level) + + # Extract country specific national MSM + PWID prevalence + iso3 <- options$area_scope + + msm_pwid_prev <- naomi.resources::load_shipp_exdata("kp_estimates", iso3) %>% + dplyr::filter(indicator == "prevalence", kp %in% c("MSM", "PWID")) + + # KP population prevalence + kp_prev <- msm_pwid_prev %>% + dplyr::select(-indicator,-lower, -upper) %>% + dplyr::mutate(median = log(median / (1 - median))) %>% + # Add in Naomi general pop prevalence + dplyr::left_join(genpop_prev, by = dplyr::join_by(area_id)) %>% + dplyr::select(kp, iso3, area_id, logit_gen_prev_msm, logit_gen_prev_pwid, median, area_level) %>% + # Calculate Log-Odds ratio + tidyr::pivot_wider(names_from = kp, + names_glue = "{.value}_{kp}", + values_from = c("median")) %>% + dplyr::mutate(msm_lor = median_MSM - logit_gen_prev_msm, + pwid_lor = median_PWID - logit_gen_prev_pwid) %>% + dplyr::select(-c("logit_gen_prev_pwid","logit_gen_prev_msm","median_PWID","median_MSM","area_level")) + + # Match KP estimates (admin0 or admin1) with SAE estimates + msm_analysis_level <- paste0("area_id",unique(msm_est$area_level)) + areas_wide <- naomi::spread_areas(areas) %>% + sf::st_drop_geometry() + map <- dplyr::select(areas_wide, area_id, dplyr::all_of(msm_analysis_level)) %>% + dplyr::rename(kp_match_area = 2) + + # Modelled estimates of proportion in each risk group + risk_group_prop <- male_srb %>% + dplyr::left_join(map, by = dplyr::join_by(area_id)) %>% + dplyr::filter(year == survey_year_sample, iso3 == options$area_scope) %>% + dplyr::select(area_id,kp_match_area, age_group, indicator, estimate_smoothed) %>% + tidyr::pivot_wider( names_from = indicator, values_from = estimate_smoothed, + values_fn = mean, names_prefix = "prop_") %>% + # Add in Naomi indicators + dplyr::right_join(naomi_est, by = c("area_id", "age_group")) %>% + dplyr::filter(!is.na(prop_nosex12m)) %>% + dplyr::mutate( + population_nosex12m = population * prop_nosex12m, + population_sexcohab = population * prop_sexcohab, + population_sexnonreg = population * prop_sexnonreg, + population_msm = population * prop_msm, + population_pwid = population * prop_pwid + ) + + # Calculate prevalence in each category + calculate_prevalence <- function(x, iso3){ + + # Log odds ratio from SRB group survey prevalence + lor <- naomi.resources:::load_shipp_exdata("srb_survey_lor", iso3) %>% + dplyr::filter(sex == "male") + + lor_15to29 <- lor$lor_15to29 + names(lor_15to29) <- lor$srb_group + + lor_30to49 <- lor$lor_30to49 + names(lor_30to49) <- lor$srb_group + + if (x$age_group[1] %in% c("Y015_019","Y020_024","Y025_029")) { + lor <- lor_15to29 + } else { + lor <- lor_30to49 + } + + population_fine <- dplyr::filter(x, indicator == "population")$estimate + plhiv <- x$plhiv[1] + kp_lor <- c("msm_lor" = x$msm_lor[1],"pwid_lor"= x$pwid_lor[1]) + lor <- c(lor[1:3], kp_lor) + prev <- logit_scale_prev(lor, population_fine, plhiv) + y <- dplyr::filter(x, indicator == "prop") %>% + dplyr::mutate(indicator = "prev",estimate = prev) + dplyr::bind_rows(x, y) + } + + + # Calculate logit prevalence and format + logit_prev <- risk_group_prop %>% + dplyr::left_join(kp_prev, by = c("kp_match_area" = "area_id")) %>% + tidyr::pivot_longer( + cols = starts_with(c("population_", "prop_")), + names_to = "indicator", + values_to = "estimate") %>% + tidyr::separate(indicator, into = c("indicator", "behav")) %>% + dplyr::filter(behav %in% c("nosex12m", "sexcohab", "sexnonreg", "msm", "pwid")) %>% + split(~ area_id + age_group) %>% + lapply(calculate_prevalence, iso3) %>% + dplyr::bind_rows() %>% + tidyr::unite("indicator", indicator, behav, sep = "_") %>% + tidyr::pivot_wider( names_from = indicator, values_from = estimate) %>% + dplyr::mutate_if(is.numeric, as.numeric) %>% + dplyr::mutate_if(is.factor, as.character) + + logit_prev + +} + + +#' Calculate the odds +#' +#' @param p Probability in [0, 1] +odds <- function(p) p / (1 - p) + +#' Calculate YWKP prevalence ratio and log odds ratio +#' +#' @param prev (General population) prevalence +#' @param fit A model relating log-odds prevalence to YWKP log odds prevalence +#' @keywords internal + +calculate_ywkp_pr_lor <- function(prev, fit = ywkp_fit) { + prev_logodds <- qlogis(prev) + prev_ywkp_logodds <- predict(fit, data.frame(prev_logodds = prev_logodds)) + # Ensure that the LOR is above that of e.g. the sexnonreg risk group + prev_ywkp_logodds <- pmax(prev_ywkp_logodds, prev_logodds + 0.25) + prev_ywkp <- plogis(prev_ywkp_logodds) + # Prevalence ratio + pr <- prev_ywkp / prev + # Log-odds ratio + lor <- prev_ywkp_logodds - prev_logodds + return(list(pr = pr, lor = lor, prev = prev, prev_ywkp = prev_ywkp)) +} + +#' Calculate prevalence and PLHIV using logit-scale disaggregation +#' +#' @param lor Log odds-ratios +#' @param N_fine Number of individuals in each group +#' @param plhiv Total number of people living with HIV +#' +#' @keywords internal + +logit_scale_prev <- function(lor, N_fine, plhiv) { + # theta represents prevalence in baseline risk group + # plogis(lor + theta) is prevalence in each risk group + # plogis(lor + theta) * N_fine is PLHIV in each risk group + optfn <- function(theta) (sum(plogis(lor + theta) * N_fine) - plhiv)^2 + # Optimisation for baseline risk group prevalence + # On the logit scale should be more numerically stable + opt <- optimise(optfn, c(-10, 10), tol = .Machine$double.eps^0.5) + # Return prevalence + plogis(lor + opt$minimum) + +} + + +#' Calculate incidence for female SRB groups. +#' +#' While maintaining age/sex/district-specific HIV incidence from Naomi, distribute +#' HIV incidence between our 4 different behavioural groups utilizing IRRs from the +#' literature +#' +#' @param naomi_output Naomi output. +#' @param options Naomi model options. +#' @param female_srb FSW adjusted estimates of female sexual risk groups generated by `shipp_adjust_sexbehav_fsw()`. +#' @param female_logit_prevalence Risk adjusted estimates of female prevalence in sexual risk groups generated by `shipp_calculate_prevalence_female()`. +#' @param survey_year Survey year to sample from the SAE model. Default is 2018. Survey year should be updated to most current household survey in the country - for countries without recent household surveys, leave at 2018 - the spatiotemporal +#' model of sexual behaviour fitted to all countries has the most data for in roughly 2018. +#' +#' @return Wide format output required for the SHIPP workbook. +#' @keywords internal + +shipp_calculate_incidence_female <- function(naomi_output, + options, + female_srb, + female_logit_prevalence, + survey_year = 2018, + kp_consensus) { + + naomi_indicators <- naomi_output %>% + dplyr::filter(indicator %in% c("population", "plhiv","prevalence","infections", "incidence"), + sex == "female", area_level == options$area_level) %>% + tidyr::pivot_wider(names_from = indicator, values_from = mean) %>% + dplyr::mutate( + incidence_cat = cut(incidence, + c(0, 0.3, 1, 3, 10^6), + labels = c("Low", "Moderate", "High", "Very High"), + include.lowest = TRUE, right = TRUE)) + + risk_group_prevalence <- female_logit_prevalence %>% + dplyr::select(area_id, age_group, starts_with("prev_")) + + df <- female_srb %>% + dplyr::filter(year == survey_year) %>% + dplyr::select(area_id, age_group, indicator, estimate_smoothed) %>% + tidyr::pivot_wider(names_from = indicator, values_from = estimate_smoothed, values_fn = mean) %>% + dplyr::left_join(naomi_indicators, by = dplyr::join_by(area_id, age_group)) %>% + dplyr::left_join(risk_group_prevalence, by = dplyr::join_by(area_id, age_group)) %>% + dplyr::filter(!is.na(population)) + + # Risk ratios for people non-regular sex partners relative to those with a + # single cohabiting sex partner + # ALPHA Network pooled analysis (Slaymaker et al CROI 2020), Jia et al systematic review, Ssempijja et al JAIDS 2022 + rr_sexcohab <- 1 + rr_sexnonreg_young <- 1.72 + rr_sexnonreg_old <- 2.1 + + # Tiered HIV risk ratio for the FSW group depending on district-level HIV + # incidence in general population + # Jones et al medRxiv "HIV incidence among women engaging in sex work in sub-Saharan Africa: a systematic review and meta-analysis" + # https://www.medrxiv.org/content/10.1101/2023.10.17.23297108v2 + # linear relationship between log(FSW incidence) and log(gen pop incidence) + # regression points shared in confidence, y = mx + b slope is 0.604104017 and + # intercept is 0.075090952 + + rr_reg_dat <- data.frame(genpop_incidence = df$incidence/100) %>% + dplyr::mutate(log_gen = log(genpop_incidence), + log_sexpaid12m = 0.604104017 * log_gen + 0.075090952, + sexpaid12m_incidence = exp(log_sexpaid12m), + rr_sexpaid12m = sexpaid12m_incidence / genpop_incidence) + + rr_sexpaid12m <- rr_reg_dat$rr_sexpaid12m + # This gives implausibly high RRs for very low districts (e.g. IRR = 297!) + # capping at 35 + rr_sexpaid12m[rr_sexpaid12m > 35] <- 35 + + # TODO: Get distributions on these and using a sampling method to get + # uncertainty in economic analysis e.g. + rr_sexnonreg_se <- 0.2 + rr_sexnonreg_se <- 1 + + + # Calculate risk group incidence + Y015_024 <- c("Y015_019", "Y020_024") + Y025_049 <- c("Y025_029","Y030_034","Y035_039","Y040_044", "Y045_049") + + df1 <- df %>% + dplyr::mutate( + rr_sexpaid12m = rr_sexpaid12m, + rr_sexnonreg = dplyr::case_when( + age_group %in% Y015_024 ~ rr_sexnonreg_young, + age_group %in% Y025_049 ~ rr_sexnonreg_old, + TRUE ~ NA_real_), + population_nosex12m = population * nosex12m, + population_sexcohab = population * sexcohab, + population_sexnonreg = population * sexnonreg, + population_sexpaid12m = population * sexpaid12m, + plhiv_nosex12m = population_nosex12m * prev_nosex12m, + plhiv_sexcohab = population_sexcohab * prev_sexcohab, + plhiv_sexnonreg = population_sexnonreg * prev_sexnonreg, + plhiv_sexpaid12m = population_sexpaid12m * prev_sexpaid12m, + susceptible_nosex12m = population_nosex12m - plhiv_nosex12m, + susceptible_sexcohab = population_sexcohab - plhiv_sexcohab, + susceptible_sexnonreg = population_sexnonreg - plhiv_sexnonreg, + susceptible_sexpaid12m = population_sexpaid12m - plhiv_sexpaid12m, + incidence_sexpaid12m = (incidence/100) * rr_sexpaid12m, + infections_sexpaid12m = susceptible_sexpaid12m * incidence_sexpaid12m) + + # Scale FSW new infections consensus estimate for KP Workbook or Goals + # Check for consensus estimate of FSW new infections + fsw_consensus <- kp_consensus[kp_consensus$key_population == "FSW", ]$infections + + if(is.na(fsw_consensus) || fsw_consensus != 0){ + # If no KP workbook present, read in FSW new infections from GOALS + goals <- naomi.resources::load_shipp_exdata("goals", "SSA") + fsw_consensus <- goals[goals$iso3 == options$area_scope, ]$`fsw-new_inf` + + } + + # Sum prior count of new infections + fsw_sum <- sum(df1$infections_sexpaid12m) + # Generate a ratio to scale FSW new infections by + fsw_ratio <- fsw_consensus / fsw_sum + + # Adjust new infections + df2 <- df1 %>% + dplyr::mutate( + # Adjust district-level new infections and incidence for FSW + infections_sexpaid12m = infections_sexpaid12m * fsw_ratio, + incidence_sexpaid12m = infections_sexpaid12m / susceptible_sexpaid12m, + # Calculate incidence in rest of groups + incidence_sexcohab = (infections - infections_sexpaid12m) / (susceptible_sexcohab + rr_sexnonreg * susceptible_sexnonreg), + incidence_sexnonreg = incidence_sexcohab * rr_sexnonreg, + infections_sexcohab = susceptible_sexcohab * incidence_sexcohab, + infections_sexnonreg = susceptible_sexnonreg * incidence_sexnonreg, + incidence_nosex12m = 0, + infections_nosex12m = 0, + rr_sexpaid12m = incidence_sexpaid12m / incidence_sexcohab + ) + + # Error here to catch that the KP adjustment has made the number of new infections + # in KPs greater than the estimated population susceptible + if(sum(df2$incidence_sexpaid12m > 1) > 0) { + stop("KP new infections exceeds susceptible population size. Please contact support.") + } + + # Calculate risk group incidence for aggregate age groups + summarise_age_cat_female <- function(dat, age_cat) { + + if (age_cat == "Y015_024") {age_groups <- c("Y015_019", "Y020_024")} + if (age_cat == "Y025_049") {age_groups <- c("Y025_029","Y030_034","Y035_039", + "Y040_044", "Y045_049")} + if (age_cat == "Y015_049") {age_groups <- c("Y015_019", "Y020_024","Y025_029", + "Y030_034","Y035_039","Y040_044", + "Y045_049")} + + + x <- dat %>% + dplyr::group_by(area_id, area_name, area_level, sex, calendar_quarter) %>% + dplyr::summarise( + "population" = sum(population * as.integer(age_group %in% age_groups)), + "plhiv" = sum(plhiv * as.integer(age_group %in% age_groups)), + "infections" = sum(infections * as.integer(age_group %in% age_groups)), + "population_nosex12m" = sum(population_nosex12m * as.integer(age_group %in% age_groups)), + "population_sexcohab" = sum(population_sexcohab * as.integer(age_group %in% age_groups)), + "population_sexnonreg" = sum(population_sexnonreg * as.integer(age_group %in% age_groups)), + "population_sexpaid12m" = sum(population_sexpaid12m * as.integer(age_group %in% age_groups)), + "plhiv_nosex12m" = sum(plhiv_nosex12m * as.integer(age_group %in% age_groups)), + "plhiv_sexnonreg" = sum(plhiv_sexnonreg * as.integer(age_group %in% age_groups)), + "plhiv_sexpaid12m" = sum(plhiv_sexpaid12m * as.integer(age_group %in% age_groups)), + "plhiv_sexcohab" = sum(plhiv_sexcohab * as.integer(age_group %in% age_groups)), + "susceptible_nosex12m" = sum(susceptible_nosex12m * as.integer(age_group %in% age_groups)), + "susceptible_sexcohab" = sum(susceptible_sexcohab * as.integer(age_group %in% age_groups)), + "susceptible_sexnonreg" = sum(susceptible_sexnonreg * as.integer(age_group %in% age_groups)), + "susceptible_sexpaid12m" = sum(susceptible_sexpaid12m * as.integer(age_group %in% age_groups)), + "infections_nosex12m" = sum(infections_nosex12m * as.integer(age_group %in% age_groups)), + "infections_sexcohab" = sum(infections_sexcohab * as.integer(age_group %in% age_groups)), + "infections_sexnonreg" = sum(infections_sexnonreg * as.integer(age_group %in% age_groups)), + "infections_sexpaid12m" = sum(infections_sexpaid12m * as.integer(age_group %in% age_groups)), + "sexnonregplus" = sum(susceptible_sexnonreg, susceptible_sexpaid12m)/(population - plhiv), + .groups = "drop") %>% + dplyr::mutate(age_group = age_cat, + nosex12m = susceptible_nosex12m/(population - plhiv), + sexcohab = susceptible_sexcohab/(population - plhiv), + sexnonreg = susceptible_sexnonreg/(population - plhiv), + sexpaid12m = susceptible_sexpaid12m/(population - plhiv), + incidence = (infections/(population - plhiv))*100, + incidence_nosex12m = infections_nosex12m/susceptible_nosex12m, + incidence_sexcohab = infections_sexcohab/susceptible_sexcohab, + incidence_sexnonreg = infections_sexnonreg/susceptible_sexnonreg, + incidence_sexpaid12m = infections_sexpaid12m/susceptible_sexpaid12m, + prev_nosex12m = plhiv_nosex12m/(susceptible_nosex12m + plhiv_nosex12m), + prev_sexcohab = plhiv_sexcohab/(susceptible_sexcohab + plhiv_sexcohab), + prev_sexnonreg = plhiv_sexnonreg/(susceptible_sexnonreg + plhiv_sexnonreg), + prev_sexpaid12m = plhiv_sexpaid12m/(susceptible_sexpaid12m + plhiv_sexpaid12m), + rr_sexpaid12m = NA) + } + + # Aggregate data + df3 <- dplyr::bind_rows(summarise_age_cat_female(df2, "Y015_024"), + summarise_age_cat_female(df2, "Y025_049"), + summarise_age_cat_female(df2, "Y015_049")) + + # Calculate incidence + df4 <- dplyr::bind_rows(df2, df3) %>% + dplyr::mutate(incidence_cat = cut(incidence, + c(0, 0.3, 1, 3, 10^6), + labels = c("Low", "Moderate", "High", "Very High"), + include.lowest = TRUE, right = TRUE)) + + # Check that sum of disaggregated infections is the same as total infections + sum_infections <- df4$infections_nosex12m + df4$infections_sexcohab + df4$infections_sexnonreg + df4$infections_sexpaid12m + + if(max(df4$infections - sum_infections) > 10^{-9}){ + stop("Risk group proportions do not sum correctly. Please contact suppport.") + } + + # Check that new infections are never negative in any behavioural risk group + sexcohab_inf_check <- sum(df4$infections_sexcohab < 0) + sexnonreg_inf_check <- sum(df4$infections_sexnonreg < 0) + sexpaid12m_inf_check <- sum(df4$infections_sexpaid12m < 0) + + if(sum(sexcohab_inf_check,sexnonreg_inf_check,sexpaid12m_inf_check)>0) { + stop("Number of new infections below 0. Please contact support.") + } + + df4 %>% + dplyr::mutate(concat = paste0(area_id, age_group), iso3 = options$area_scope) %>% + dplyr::select(area_id, age_group, concat, + nosex12m, sexcohab, sexnonregplus, sexnonreg, sexpaid12m, + iso3, area_level, + population, plhiv, infections, incidence, incidence_cat, + prev_nosex12m, prev_sexcohab, prev_sexnonreg, prev_sexpaid12m, + rr_sexpaid12m, # rr_sexnonreg, + population_nosex12m, population_sexcohab, + population_sexnonreg, population_sexpaid12m, + plhiv_nosex12m, plhiv_sexcohab, + plhiv_sexnonreg, plhiv_sexpaid12m, + susceptible_nosex12m, susceptible_sexcohab, + susceptible_sexnonreg, susceptible_sexpaid12m, + incidence_nosex12m, incidence_sexcohab, + incidence_sexnonreg, incidence_sexpaid12m, + infections_nosex12m,infections_sexcohab, + infections_sexnonreg, infections_sexpaid12m, + incidence_cat) %>% + dplyr::mutate_if(is.numeric, as.numeric) %>% + dplyr::mutate_if(is.factor, as.character) + +} + +#' Calculate incidence in all behavioural groups +#' +#' @param outputs Naomi output. +#' @param options Naomi model options. +#' @param male_srb MSM and PWID adjusted estimated of male SRB groups generated by `shipp_adjust_sexbehav_msm_pwid()`. +#' @param male_logit_prevalence Risk adjusted estimates of male HIV prevalence in sexual risk groups generated by `shipp_calculate_prevalence_male()`. +#' @param survey_year Survey year to sample from the SAE model. Default is 2018. Survey year should be updated to most current household survey in the country - for countries without recent household surveys, leave at 2018 - the spatiotemporal +#' model of sexual behaviour fitted to all countries has the most data for in roughly 2018. + +#' @return Wide format output required for the SHIPP workbook + +shipp_calculate_incidence_male <- function(naomi_output, + options, + male_srb, + male_logit_prevalence, + survey_year = 2018, + kp_consensus) { + + + naomi_indicators <- naomi_output %>% + dplyr::filter(indicator %in% c("population", "plhiv","prevalence","infections", "incidence"), + sex == "male", area_level == options$area_level) %>% + tidyr::pivot_wider(names_from = indicator, values_from = mean) %>% + dplyr::mutate( + incidence_cat = cut(incidence, + c(0, 0.3, 1, 3, 10^6), + labels = c("Low", "Moderate", "High", "Very High"), + include.lowest = TRUE, right = TRUE)) + + risk_group_prevalence <- male_logit_prevalence %>% + dplyr::select(area_id, age_group, gen_prev, starts_with( "prev_")) %>% + dplyr::mutate(msm_pr = prev_msm/ gen_prev, + pwid_pr = prev_pwid / gen_prev) + + df <- male_srb %>% + dplyr::filter(year == survey_year) %>% + dplyr::select(area_id, age_group, indicator, estimate_smoothed) %>% + tidyr::pivot_wider(names_from = indicator, values_from = estimate_smoothed, values_fn = mean) %>% + dplyr::left_join(naomi_indicators, by = dplyr::join_by(area_id, age_group)) %>% + dplyr::left_join(risk_group_prevalence, by = dplyr::join_by(area_id, age_group)) %>% + dplyr::filter(!is.na(population)) + + + # ALPHA Network pooled analysis (Slaymaker et al CROI 2020), Hoffman et al JAIDS 2022, Ssempijja et al JAIDS 2022 + rr_sexcohab <- 1 + rr_sexnonreg_young <- 1.89 + rr_sexnonreg_old <- 2.1 + + # TODO: Get distributions on these and using a sampling method to get uncertainty in economic analysis e.g. + rr_sexnonreg_se <- 0.2 + rr_sexnonreg_se <- 1 + + + # Calculate risk group incidence + Y015_024 <- c("Y015_019", "Y020_024") + Y025_049 <- c("Y025_029","Y030_034","Y035_039","Y040_044", "Y045_049") + + df1 <- df %>% + dplyr::mutate( + msm_pr = round(prev_msm / prevalence, 2), + pwid_pr = round(prev_pwid / prevalence, 2), + # Setting artificial cutoff of IRR of 2.5 due to Stannah et al + # Lancet HIV systematic review of MSM vs gen pop IRR + # https://www.thelancet.com/journals/lanhiv/article/PIIS2352-3018(23)00111-X/fulltext + rr_msm = dplyr::if_else(msm_pr > 2.5, msm_pr, 2.5), + rr_pwid = dplyr::if_else(pwid_pr > 2.5, pwid_pr, 2.5), + rr_sexnonreg = dplyr::case_when( + age_group %in% Y015_024 ~ rr_sexnonreg_young, + age_group %in% Y025_049 ~ rr_sexnonreg_old, + TRUE ~ NA_real_), + population_nosex12m = population * nosex12m, + population_sexcohab = population * sexcohab, + population_sexnonreg = population * sexnonreg, + population_msm = population * msm, + population_pwid = population * pwid, + plhiv_nosex12m = population_nosex12m * prev_nosex12m, + plhiv_sexcohab = population_sexcohab * prev_sexcohab, + plhiv_sexnonreg = population_sexnonreg * prev_sexnonreg, + plhiv_msm = population_msm * prev_msm, + plhiv_pwid = population_pwid * prev_pwid, + susceptible_nosex12m = population_nosex12m - plhiv_nosex12m, + susceptible_sexcohab = population_sexcohab - plhiv_sexcohab, + susceptible_sexnonreg = population_sexnonreg - plhiv_sexnonreg, + susceptible_msm = population_msm - plhiv_msm, + susceptible_pwid = population_pwid - plhiv_pwid, + incidence_msm = (incidence/100) * rr_msm, + incidence_pwid = (incidence/100) * rr_pwid, + infections_msm = susceptible_msm * incidence_msm, + infections_pwid = susceptible_pwid * incidence_pwid) + + # Scale MSM and PWID new infections consensus estimate for KP Workbook or Goals + # Check for consensus estimate of MSM and PWID new infections + msm_consensus <- kp_consensus[kp_consensus$key_population == "MSM", ]$infections + pwid_consensus <- kp_consensus[kp_consensus$key_population == "PWID", ]$infections + + if(is.na(msm_consensus) || msm_consensus == 0){ + # If no KP workbook present, read in FSW new infections from GOALS + goals <- naomi.resources::load_shipp_exdata("goals", "SSA") + msm_consensus <- goals[goals$iso3 == options$area_scope, ]$`msm-new_inf` + } + + if(is.na(pwid_consensus) || pwid_consensus == 0){ + # If no KP workbook present, read in FSW new infections from GOALS + goals <- naomi.resources::load_shipp_exdata("goals", "SSA") + pwid_consensus <- goals[goals$iso3 == options$area_scope, ]$`pwid-new_inf` + } + + # scale consensus that we'll use to account for the 1:10 ratio assumption of + # male:female PWID + pwid_consensus <- pwid_consensus * 0.91 + + # Sum prior count of new infections + msm_sum <- sum(df1$infections_msm) + pwid_sum <- sum(df1$infections_pwid) + + # Generate a ratio to scale MSM and PWID new infections by + msm_ratio <- msm_consensus / msm_sum + pwid_ratio <- pwid_consensus / pwid_sum + + # Adjust new infections + df2 <- df1 %>% + dplyr::mutate( + # Adjust district-level new infections and incidence for MSM + infections_msm = infections_msm * msm_ratio, + incidence_msm = infections_msm / susceptible_msm, + # Adjust district-level new infections and incidence for PWID + infections_pwid = infections_pwid * pwid_ratio, + incidence_pwid = infections_pwid / susceptible_pwid, + + # Adjust sexcohab and sexnonreg new infections and incidence to scale rest of infections + # from the district + incidence_nosex12m = 0, + incidence_sexcohab = (infections - infections_msm - infections_pwid) / (susceptible_sexcohab + + rr_sexnonreg * susceptible_sexnonreg), + incidence_sexnonreg = incidence_sexcohab * rr_sexnonreg, + infections_nosex12m = 0, + infections_sexcohab = susceptible_sexcohab * incidence_sexcohab, + infections_sexnonreg = susceptible_sexnonreg * incidence_sexnonreg, + + rr_msm = incidence_msm / incidence_sexcohab, + rr_pwid = incidence_pwid / incidence_sexcohab + ) + + # Error here to catch that the KP adjustment has made the number of new infections + # in KPs greater than the estimated population susceptible + if(sum(df2$incidence_msm > 1, df2$incidence_pwid > 1) > 0) { + stop("KP new infections exceeds susceptible population size. Please contact support.") + } + + + # Calculate risk group incidence for aggregate age groups + + summarise_age_cat_male <- function(dat, age_cat) { + + if (age_cat == "Y015_024") {age_groups <- c("Y015_019", "Y020_024")} + if (age_cat == "Y025_049") {age_groups <- c("Y025_029","Y030_034","Y035_039", + "Y040_044", "Y045_049")} + if (age_cat == "Y015_049") {age_groups <- c("Y015_019", "Y020_024","Y025_029", + "Y030_034","Y035_039","Y040_044", + "Y045_049")} + + dat %>% + dplyr::group_by(area_id, area_name, area_level, sex, calendar_quarter) %>% + dplyr::summarise( + "population" = sum(population * as.integer(age_group %in% age_groups)), + "plhiv" = sum(plhiv * as.integer(age_group %in% age_groups)), + "infections" = sum(infections * as.integer(age_group %in% age_groups)), + "population_nosex12m" = sum(population_nosex12m * as.integer(age_group %in% age_groups)), + "population_sexcohab" = sum(population_sexcohab * as.integer(age_group %in% age_groups)), + "population_sexnonreg" = sum(population_sexnonreg * as.integer(age_group %in% age_groups)), + "population_msm" = sum(population_pwid * as.integer(age_group %in% age_groups)), + "population_pwid" = sum(population_pwid * as.integer(age_group %in% age_groups)), + "plhiv_nosex12m" = sum(plhiv_nosex12m * as.integer(age_group %in% age_groups)), + "plhiv_sexnonreg" = sum(plhiv_sexnonreg * as.integer(age_group %in% age_groups)), + "plhiv_sexcohab" = sum(plhiv_sexcohab * as.integer(age_group %in% age_groups)), + "plhiv_msm" = sum(plhiv_msm * as.integer(age_group %in% age_groups)), + "plhiv_pwid" = sum(plhiv_pwid * as.integer(age_group %in% age_groups)), + "susceptible_nosex12m" = sum(susceptible_nosex12m * as.integer(age_group %in% age_groups)), + "susceptible_sexcohab" = sum(susceptible_sexcohab * as.integer(age_group %in% age_groups)), + "susceptible_sexnonreg" = sum(susceptible_sexnonreg * as.integer(age_group %in% age_groups)), + "susceptible_msm" = sum(susceptible_msm * as.integer(age_group %in% age_groups)), + "susceptible_pwid" = sum(susceptible_pwid * as.integer(age_group %in% age_groups)), + "infections_nosex12m" = sum(infections_nosex12m * as.integer(age_group %in% age_groups)), + "infections_sexcohab" = sum(infections_sexcohab * as.integer(age_group %in% age_groups)), + "infections_sexnonreg" = sum(infections_sexnonreg * as.integer(age_group %in% age_groups)), + "infections_msm" = sum(infections_msm * as.integer(age_group %in% age_groups)), + "infections_pwid" = sum(infections_pwid * as.integer(age_group %in% age_groups)), + "sexnonregplus" = sum(susceptible_sexnonreg)/(population - plhiv), + .groups = "drop") %>% + dplyr::mutate(age_group = age_cat, + nosex12m = susceptible_nosex12m/(population - plhiv), + sexcohab = susceptible_sexcohab/(population - plhiv), + sexnonreg = susceptible_sexnonreg/(population - plhiv), + msm = susceptible_msm/(population - plhiv), + pwid = susceptible_pwid/(population - plhiv), + incidence = (infections/(population - plhiv)) * 100, + incidence_nosex12m = infections_nosex12m/susceptible_nosex12m, + incidence_sexcohab = infections_sexcohab/susceptible_sexcohab, + incidence_sexnonreg = infections_sexnonreg/susceptible_sexnonreg, + incidence_msm = infections_msm/susceptible_msm, + incidence_pwid = infections_pwid/susceptible_pwid, + prev_nosex12m = plhiv_nosex12m/(susceptible_nosex12m + plhiv_nosex12m), + prev_sexcohab = plhiv_sexcohab/(susceptible_sexcohab + plhiv_sexcohab), + prev_sexnonreg = plhiv_sexnonreg/(susceptible_sexnonreg + plhiv_sexnonreg), + prev_msm = plhiv_msm/(susceptible_msm + plhiv_msm), + prev_pwid = plhiv_pwid/(susceptible_pwid + plhiv_pwid), + rr_msm = NA, + rr_pwid = NA) + } + + # Aggregate data + df3 <- dplyr::bind_rows(summarise_age_cat_male(df2, "Y015_024"), + summarise_age_cat_male(df2, "Y025_049"), + summarise_age_cat_male(df2, "Y015_049")) + + # Calculate incidence + df4 <- dplyr::bind_rows(df2, df3) %>% + dplyr::mutate(incidence_cat = cut(incidence, + c(0, 0.3, 1, 3, 10^6), + labels = c("Low", "Moderate", "High", "Very High"), + include.lowest = TRUE, right = TRUE)) + + + + # Check that sum of disaggregated infections is the same as total infections + # TO DO: add warning for sum not matching - contact admin + sum_infections <- df4$infections_nosex12m + df4$infections_sexcohab + df4$infections_sexnonreg + df4$infections_msm + df4$infections_pwid + + if(max(df4$infections - sum_infections) > 10^{-9}){ + stop("Risk group proportions do not sum correctly. Please contact suppport.") + } + + # Check that new infections are never negative in any behavioural risk group + sexcohab_inf_check <- sum(df4$infections_sexcohab < 0) + sexnonreg_inf_check <- sum(df4$infections_sexnonreg < 0) + msm_inf_check <- sum(df4$infections_msm < 0) + pwid_inf_check <- sum(df4$infections_pwid < 0) + if(sum(sexcohab_inf_check,sexnonreg_inf_check,msm_inf_check,pwid_inf_check)>0) { + stop("Number of new infections below 0. Please contact support.") + } + + + df4 %>% + dplyr::mutate(concat = paste0(area_id, age_group), iso3 = options$area_scope) %>% + dplyr::select(area_id, age_group, concat, + nosex12m, sexcohab, sexnonregplus, sexnonreg, msm, pwid, + iso3, area_level, + population, plhiv, infections, incidence, incidence_cat, + prev_nosex12m, prev_sexcohab, + prev_sexnonreg, prev_msm, prev_pwid, + rr_msm, rr_pwid, rr_sexnonreg, + population_nosex12m, population_sexcohab, + population_sexnonreg, population_msm, population_pwid, + plhiv_nosex12m, plhiv_sexcohab, + plhiv_sexnonreg, plhiv_msm, plhiv_pwid, + susceptible_nosex12m, susceptible_sexcohab, + susceptible_sexnonreg, susceptible_msm, susceptible_pwid, + incidence_nosex12m, incidence_sexcohab, + incidence_sexnonreg, incidence_msm, incidence_pwid, + infections_nosex12m,infections_sexcohab, + infections_sexnonreg, infections_msm, infections_pwid) %>% + dplyr::mutate_if(is.numeric, as.numeric) %>% + dplyr::mutate_if(is.factor, as.character) + + + +} + + +#' Generate outputs to update SHIPP tool. +#' +#' @param naomi_output Path to naomi output (zip file or hintr object). +#' @param pjnz Path to spectrum file. +#' @param male_srb Estimates of male sexual risk groups generated by `shipp_adjust_sexbehav_msm_pwid()` +#' @param male_logit_prevalence Risk adjusted estimates of male prevalence in sexual risk groups generated by `shipp_calculate_prevalence_male()` +#' @param survey_year Survey year to sample from the SAE model. Default is 2018. Survey year should be updated to most current household survey in the country - for countries without recent household surveys, leave at 2018 - the spatiotemporal +#' model of sexual behaviour fitted to all countries has the most data for in roughly 2018. +#' +#' @return Output files to update SHIPP excel workbook. +#' @keywords internal + + +shipp_generate_risk_populations <- function(naomi_output, + pjnz, + survey_year = 2018) { + + # Read in naomi outputs + if (tolower(tools::file_ext(naomi_output)) %in% c("rds", "qs")) { + # Read files if hintr rds provided + model_object <- read_hintr_output(naomi_output) + outputs <- model_object$output_package + options <- outputs$fit$model_options + + } else if (grepl("\\.zip$", naomi_output)) { + # Read files if output zip is provided + output_zip <- naomi_output + outputs <- naomi::read_output_package(output_zip) + options <- outputs$fit$model_options + } + + # Check for concordence between area_ids in shipp resources from `naomi.resources` + # and Naomi estimates + assert_shipp_resource_hierarchy(outputs, options) + + # Format naomi output + naomi <- shipp_format_naomi(outputs, options) + + # Naomi population + naomi_pop <- naomi$naomi_long %>% + dplyr::filter(indicator == "population") %>% + dplyr::select(area_id, area_level,sex, age_group, area_level, + spectrum_region_code, population = mean) + + naomi_pop$iso3 <- options$area_scope + + # Disaggregate KP PSEs from Oli's analysis to 5-year bands + kp_consensus <- extract_kp_workbook(pjnz) + fsw_est <- shipp_disaggregate_fsw(outputs, options, naomi_pop, kp_consensus) + pwid_est <- shipp_disaggregate_pwid(outputs, options, naomi_pop, kp_consensus) + msm_est <- shipp_disaggregate_msm(outputs, options, naomi_pop, kp_consensus) + + # Adjust SAE model output with KP proportions + female_srb <- shipp_adjust_sexbehav_fsw(outputs, options, fsw_est) + male_srb <- shipp_adjust_sexbehav_msm_pwid(outputs, options, msm_est, pwid_est) + + # Calculate risk group prevalence + female_logit_prevalence <- shipp_calculate_prevalence_female(naomi$naomi_long, + options, + fsw_est, + female_srb, + survey_year) + + male_logit_prevalence <- shipp_calculate_prevalence_male(naomi$naomi_long, + outputs$meta_area, + options, + msm_est, + male_srb, + survey_year) + + # Calculate risk group incidence + female_incidence <- shipp_calculate_incidence_female(naomi$naomi_long, + options, + female_srb, + female_logit_prevalence, + survey_year, + kp_consensus) + + male_incidence <- shipp_calculate_incidence_male(naomi$naomi_long, + options, + male_srb, + male_logit_prevalence, + survey_year, + kp_consensus) + + meta <- data.frame(kp = c("FSW", "MSM", "PWID"), + consensus_estimate = c(unique(fsw_est$consensus_estimate), + unique(msm_est$consensus_estimate), + unique(pwid_est$consensus_estimate))) + + + + v <- list(female_incidence = female_incidence, + male_incidence = male_incidence, + naomi_output = naomi$naomi_wide, + meta_consensus = meta) + + v + +} + +#' Throw warning when area hierarchy in external SHIPP resources read in from +#' `naomi.resources` do not match Naomi outputs used to update SHIPP estimates. +#' +#' @param outputs Naomi outputs. +#' @param options Naomi options. +#' @return Resource hierarchy +#' @keywords internal + + +assert_shipp_resource_hierarchy <- function(outputs, + options){ + + # iso3 from model options + iso3 <- options$area_scope + + # KP PSE's + pse <- naomi.resources::load_shipp_exdata("kp_estimates", iso3) + + # SRB SAE model estimates + female_srb <- naomi.resources::load_shipp_exdata("srb_female", options$area_scope) + male_srb <- naomi.resources::load_shipp_exdata("srb_male", options$area_scope) + + # Naomi area_ids at lowest admin-level + naomi_hierarchy <- outputs$meta_area %>% dplyr::filter(area_level == options$area_level) + naomi_ids <- unique(naomi_hierarchy$area_id) + + pse_diff <- setdiff(unique(pse$area_id), unique(outputs$meta_area$area_id)) + female_srb_diff <- setdiff(unique(female_srb$area_id), naomi_ids) + male_srb_diff <- setdiff(unique(male_srb$area_id), naomi_ids) + + if(length(pse_diff) != 0 ){ + + stop(paste0(t_("SHIPP_ERROR_KP_PREFIX"), + paste0(unique(pse$area_id), collapse = "; "), + t_("SHIPP_ERROR_NAOMI_MISMATCH"), + paste0(unique(naomi_ids), collapse = "; "), + t_("SHIPP_ERROR_CONTACT_SUPPORT"))) + } + + if(length(female_srb_diff) != 0 ){ + + stop(paste0(t_("SHIPP_ERROR_FSRB_PREFIX"), + paste0(unique(pse$area_id), collapse = "; "), + t_("SHIPP_ERROR_NAOMI_MISMATCH"), + paste0(unique(naomi_ids), collapse = "; "), + t_("SHIPP_ERROR_CONTACT_SUPPORT"))) + } + + if(length(male_srb_diff) != 0 ){ + + stop(paste0(t_("SHIPP_ERROR_MSRB_PREFIX"), + paste0(unique(pse$area_id), collapse = "; "), + t_("SHIPP_ERROR_NAOMI_MISMATCH"), + paste0(unique(naomi_ids), collapse = "; "), + t_("SHIPP_ERROR_CONTACT_SUPPORT"))) + } + +} + +write_shipp_workbook <- function(sheets, dest) { + template_path <- naomi.resources::get_shipp_workbook_path() + withCallingHandlers( + write_xlsx_sheets(template_path, sheets, dest), + error = function(e) { + e$message <- t_("SHIPP_ERROR_WRITE", list(message = e$message)) + stop(e) + }) +} diff --git a/R/test-helpers.R b/R/test-helpers.R new file mode 100644 index 00000000..acd6e7a9 --- /dev/null +++ b/R/test-helpers.R @@ -0,0 +1,49 @@ +## File contains test helpers which we want to use both here and in hintr + +#' Build JSON from template and a set of params +#' +#' @param naomi_output Calibrated naomi output +#' +#' @return Calibrated naomi output matched to MWI test data on +#' `naomi.resources` to be used to generate the shipp tool. +#' @keywords internal +make_shipp_testfiles <- function(naomi_output) { + # Create naomi outputs align with testing data in naomi.resources: + # - Change iso3 to "MWI_demo" + # - Restrict outputs to admin2 + output <- naomi::read_hintr_output(naomi_output$model_output_path) + + # Areas + meta_area_demo <- dplyr::mutate(output$output_package$meta_area, + area_id = dplyr::if_else(area_id == "MWI", "MWI_demo", area_id), + parent_area_id = dplyr::if_else(parent_area_id == "MWI", "MWI_demo", parent_area_id)) + + meta_area_demo <- dplyr::filter(meta_area_demo, area_level <= 2) + + # Indicators + ind_demo <- dplyr::mutate(output$output_package$indicators, + area_id = dplyr::if_else(area_id == "MWI", "MWI_demo", area_id)) + + ind_demo <- dplyr::filter(ind_demo, area_id %in% meta_area_demo$area_id) + + + # Options + options_demo <- output$output_package$fit$model_options + options_demo$area_scope <- "MWI_demo" + options_demo$area_level <- 2 + + # Save out demo output package + demo <- output + demo$output_package$indicators <- ind_demo + demo$output_package$fit$model_options <- options_demo + demo$output_package$meta_area <- meta_area_demo + + out_demo <- tempfile(fileext = ".qs") + naomi:::hintr_save(demo, out_demo) + + # Add to existing hintr_test data + shipp_output_demo <- naomi_output + shipp_output_demo$model_output_path <- out_demo + + shipp_output_demo +} diff --git a/R/utils.R b/R/utils.R index 6662f84c..94f2f598 100644 --- a/R/utils.R +++ b/R/utils.R @@ -97,6 +97,38 @@ vlapply <- function(X, FUN, ...) { vapply(X, FUN, ..., FUN.VALUE = logical(1)) } +vnapply <- function(X, FUN, ...) { + vapply(X, FUN, ..., FUN.VALUE = numeric(1)) +} + +vcapply <- function(X, FUN, ...) { + vapply(X, FUN, ..., FUN.VALUE = character(1)) +} + is_empty <- function(x) { length(x) == 0 || is.null(x) || is.na(x) || !nzchar(x) } + +#' Write list of data frames into an xlsx file +#' +#' This doesn't write colmn headers into the workbook, it expects +#' that these already exist. +#' +#' @param template Path to xlsx file with sheets +#' @param sheets Named list of data frames to write into template. The names +#' must match the destination sheet in the xlsx +#' @param path Path to output the filled in xlsx +#' +#' @return Path to complete xlsx file +#' @keywords internal +write_xlsx_sheets <- function(template, sheets, path) { + wb <- openxlsx::loadWorkbook(template) + for (sheet in names(sheets)) { + openxlsx::writeData(wb, sheet, sheets[[sheet]], + startRow = 2, + colNames = FALSE) + } + + openxlsx::saveWorkbook(wb, path) + path +} diff --git a/docker/Dockerfile b/docker/Dockerfile index 99299213..f5a05a6e 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -13,7 +13,11 @@ RUN install_packages --repo=https://mrc-ide.r-universe.dev \ mockr \ rvest \ pkgbuild \ - testthat.buildkite + testthat.buildkite \ + openxlsx + +RUN install_remote \ + mrc-ide/naomi.resources ## Model run will try to parallelise over as many threads as are available ## potentially slowing the application, manually limit threads to 1 diff --git a/inst/traduire/en-translation.json b/inst/traduire/en-translation.json index e33afc4b..6ed50714 100644 --- a/inst/traduire/en-translation.json +++ b/inst/traduire/en-translation.json @@ -91,7 +91,7 @@ "PROGRESS_DOWNLOAD_COARSE": "Generating coarse-output download", "PROGRESS_DOWNLOAD_SUMMARY": "Generating summary report", "PROGRESS_DOWNLOAD_COMPARISON": "Generating comparison report", - "PROGRESS_DOWNLOAD_AGYW": "Generating AGYW tool", + "PROGRESS_DOWNLOAD_SHIPP": "Generating SHIPP tool", "PLHIV_MAP_PLOT_TITLE": "People living with HIV (15+)", "ART_CURRENT_MAP_PLOT_TITLE": "Residents receiving ART (15+)", "INFECTIONS_MAP_PLOT_TITLE": "Annual HIV infections (15+)", @@ -275,5 +275,14 @@ "DOWNLOAD_OUTPUT_DESCRIPTION": "Naomi output uploaded from Naomi web app", "DOWNLOAD_SUMMARY_DESCRIPTION": "Naomi summary report uploaded from Naomi web app", "DOWNLOAD_COMPARISON_DESCRIPTION": "Naomi comparison report uploaded from Naomi web app", - "DOWNLOAD_AGYW_DESCRIPTION": "Naomi AGYW tool uploaded from Naomi web app" + "DOWNLOAD_SHIPP_DESCRIPTION": "Naomi SHIPP tool uploaded from Naomi web app", + "SHIPP_ERROR_KP_PREFIX": "Available KP PSE estimates for: \n", + "SHIPP_ERROR_FSRB_PREFIX": " Available female sexual risk behaviour survey estimates for: \n", + "SHIPP_ERROR_MSRB_PREFIX": " Available male sexual risk behaviour survey estimates for: \n", + "SHIPP_ERROR_NAOMI_MISMATCH": "\n\n Do not match Naomi estimates for: \n", + "SHIPP_ERROR_CONTACT_SUPPORT": "\n\nTo update estimates, please contact Naomi support.", + "SHIPP_ERROR_WRITE": "Failed to build workbook, please contact support: {{message}}", + "FSW_CONSENSUS_WARNING": "Consensus estimate for FSW is larger than 5% of total females aged 15-59. Modelled estimates for FWS PSEs will be used.", + "MSM_CONSENSUS_WARNING": "Consensus estimate for MSM is larger than 5% of total males aged 15-59. Modelled estimates for FWS PSEs will be used.", + "PWID_CONSENSUS_WARNING": "Consensus estimate for PWID is larger than 5% of total males aged 15-59. Modelled estimates for PWID PSEs will be used." } diff --git a/inst/traduire/fr-translation.json b/inst/traduire/fr-translation.json index 7bc5a32b..ae34a03f 100644 --- a/inst/traduire/fr-translation.json +++ b/inst/traduire/fr-translation.json @@ -139,6 +139,7 @@ "PROGRESS_DOWNLOAD_COARSE": "Génération d'un téléchargement à sortie grossière", "PROGRESS_DOWNLOAD_SUMMARY": "Générer un rapport de synthèse", "PROGRESS_DOWNLOAD_COMPARISON": "Générer un rapport de comparaison", + "PROGRESS_DOWNLOAD_SHIPP": "Generating SHIPP tool", "CANNOT_RECALIBRATE": "L'étalonnage ne peut pas être relancé pour cet ajustement de modèle, veuillez relancer l'étape d'ajustement.", "INVALID_ART_TIME_PERIOD": "Données trimestrielles non fournies pour tous les désagrégats ou données annuelles dupliquées fournies pour tous les désagrégats.", "ART_TOTAL": "Numéros sur TARV", @@ -273,5 +274,15 @@ "BIRTHS_FACILITY_DESC": "Nombre de naissances dans les établissements de santé", "DOWNLOAD_OUTPUT_DESCRIPTION": "Paquet Naomi téléchargée depuis l'application web Naomi", "DOWNLOAD_SUMMARY_DESCRIPTION": "Rapport de synthèse Naomi téléchargé depuis l'application web Naomi", - "DOWNLOAD_COMPARISON_DESCRIPTION": "Rapport de comparaison Naomi téléchargé à partir de l'application web Naomi" + "DOWNLOAD_COMPARISON_DESCRIPTION": "Rapport de comparaison Naomi téléchargé à partir de l'application web Naomi", + "DOWNLOAD_SHIPP_DESCRIPTION": "Naomi SHIPP tool uploaded from Naomi web app", + "SHIPP_ERROR_KP_PREFIX": " Available KP PSE estimates for: \n", + "SHIPP_ERROR_FSRB_PREFIX": " Available female sexual risk behaviour survey estimates for: \n", + "SHIPP_ERROR_MSRB_PREFIX": " Available male sexual risk behaviour survey estimates for: \n", + "SHIPP_ERROR_NAOMI_MISMATCH": "\n\n Do not match Naomi estimates for: \n", + "SHIPP_ERROR_CONTACT_SUPPORT": "\n\nTo update estimates, please contact Naomi support.", + "SHIPP_ERROR_WRITE": "Échec de la création du classeur, veuillez contacter l'assistance: {{message}}", + "FSW_CONSENSUS_WARNING": "Consensus estimate for FSW is larger than 5% of total females aged 15-59. Modelled estimates for FWS PSEs will be used.", + "MSM_CONSENSUS_WARNING": "Consensus estimate for MSM is larger than 5% of total males aged 15-59. Modelled estimates for FWS PSEs will be used.", + "PWID_CONSENSUS_WARNING": "Consensus estimate for PWID is larger than 5% of total males aged 15-59. Modelled estimates for PWID PSEs will be used." } diff --git a/inst/traduire/pt-translation.json b/inst/traduire/pt-translation.json index 07ed65c2..c23747e5 100644 --- a/inst/traduire/pt-translation.json +++ b/inst/traduire/pt-translation.json @@ -139,6 +139,7 @@ "PROGRESS_DOWNLOAD_COARSE": "Geração de resultados grosseiros", "PROGRESS_DOWNLOAD_SUMMARY": "Elaboração de relatório de síntese", "PROGRESS_DOWNLOAD_COMPARISON": "Elaboração de relatório de comparação", + "PROGRESS_DOWNLOAD_SHIPP": "Generating SHIPP tool", "CANNOT_RECALIBRATE": "A calibração não pode ser executada de novo para este modelo de ajuste, por favor, execute de novo o passo de ajuste.", "INVALID_ART_TIME_PERIOD": "Dados trimestrais não fornecidos para todos os desagregados ou dados anuais duplicados fornecidos para todos os desagregados.", "ART_TOTAL": "Contagem de TARV", @@ -273,5 +274,15 @@ "BIRTHS_FACILITY_DESC": "Número de nascimentos em estabelecimentos de saúde", "DOWNLOAD_OUTPUT_DESCRIPTION": "Pacote Naomi descarregado a partir da aplicação web Naomi", "DOWNLOAD_SUMMARY_DESCRIPTION": "Relatório de síntese da Naomi carregado da aplicação web Naomi", - "DOWNLOAD_COMPARISON_DESCRIPTION": "Relatório de comparação Naomi carregado a partir da aplicação web Naomi" + "DOWNLOAD_COMPARISON_DESCRIPTION": "Relatório de comparação Naomi carregado a partir da aplicação web Naomi", + "DOWNLOAD_SHIPP_DESCRIPTION": "Naomi SHIPP tool uploaded from Naomi web app", + "SHIPP_ERROR_KP_PREFIX": " Available KP PSE estimates for: \n", + "SHIPP_ERROR_FSRB_PREFIX": " Available female sexual risk behaviour survey estimates for: \n", + "SHIPP_ERROR_MSRB_PREFIX": " Available male sexual risk behaviour survey estimates for: \n", + "SHIPP_ERROR_NAOMI_MISMATCH": "\n\n Do not match Naomi estimates for: \n", + "SHIPP_ERROR_CONTACT_SUPPORT": "\n\nTo update estimates, please contact Naomi support.", + "SHIPP_ERROR_WRITE": "Falha ao criar a pasta de trabalho, entre em contato com o suporte: {{message}}", + "FSW_CONSENSUS_WARNING": "Consensus estimate for FSW is larger than 5% of total females aged 15-59. Modelled estimates for FWS PSEs will be used.", + "MSM_CONSENSUS_WARNING": "Consensus estimate for MSM is larger than 5% of total males aged 15-59. Modelled estimates for FWS PSEs will be used.", + "PWID_CONSENSUS_WARNING": "Consensus estimate for PWID is larger than 5% of total males aged 15-59. Modelled estimates for PWID PSEs will be used." } diff --git a/man/assert_shipp_resource_hierarchy.Rd b/man/assert_shipp_resource_hierarchy.Rd new file mode 100644 index 00000000..f4e3f133 --- /dev/null +++ b/man/assert_shipp_resource_hierarchy.Rd @@ -0,0 +1,15 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{assert_shipp_resource_hierarchy} +\alias{assert_shipp_resource_hierarchy} +\title{Throw warning when area hierarchy in external SHIPP resources read in from +\code{naomi.resources} do not match Naomi outputs used to update SHIPP estimates.} +\usage{ +assert_shipp_resource_hierarchy(outputs, options) +} +\description{ +@param outputs Naomi outputs. +@param options Naomi options. +@return Resource hierarchy +@keywords internal +} diff --git a/man/calculate_ywkp_pr_lor.Rd b/man/calculate_ywkp_pr_lor.Rd new file mode 100644 index 00000000..b481a540 --- /dev/null +++ b/man/calculate_ywkp_pr_lor.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{calculate_ywkp_pr_lor} +\alias{calculate_ywkp_pr_lor} +\title{Calculate YWKP prevalence ratio and log odds ratio} +\usage{ +calculate_ywkp_pr_lor(prev, fit = ywkp_fit) +} +\arguments{ +\item{prev}{(General population) prevalence} + +\item{fit}{A model relating log-odds prevalence to YWKP log odds prevalence} +} +\description{ +Calculate YWKP prevalence ratio and log odds ratio +} +\keyword{internal} diff --git a/man/hintr_prepare_agyw_download.Rd b/man/hintr_prepare_shipp_download.Rd similarity index 57% rename from man/hintr_prepare_agyw_download.Rd rename to man/hintr_prepare_shipp_download.Rd index ecf3048f..3d5d40fd 100644 --- a/man/hintr_prepare_agyw_download.Rd +++ b/man/hintr_prepare_shipp_download.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/downloads.R -\name{hintr_prepare_agyw_download} -\alias{hintr_prepare_agyw_download} -\title{Prepare AGYW tool download} +\name{hintr_prepare_shipp_download} +\alias{hintr_prepare_shipp_download} +\title{Prepare SHIPP tool download} \usage{ -hintr_prepare_agyw_download(output, pjnz, path = tempfile(fileext = ".xlsx")) +hintr_prepare_shipp_download(output, pjnz, path = tempfile(fileext = ".xlsx")) } \arguments{ \item{pjnz}{Path to input PJNZ file} @@ -17,5 +17,5 @@ hintr_prepare_agyw_download(output, pjnz, path = tempfile(fileext = ".xlsx")) Path to output file and metadata for file } \description{ -Prepare AGYW tool download +Prepare SHIPP tool download } diff --git a/man/logit_scale_prev.Rd b/man/logit_scale_prev.Rd new file mode 100644 index 00000000..8f3f7feb --- /dev/null +++ b/man/logit_scale_prev.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{logit_scale_prev} +\alias{logit_scale_prev} +\title{Calculate prevalence and PLHIV using logit-scale disaggregation} +\usage{ +logit_scale_prev(lor, N_fine, plhiv) +} +\arguments{ +\item{lor}{Log odds-ratios} + +\item{N_fine}{Number of individuals in each group} + +\item{plhiv}{Total number of people living with HIV} +} +\description{ +Calculate prevalence and PLHIV using logit-scale disaggregation +} +\keyword{internal} diff --git a/man/make_shipp_testfiles.Rd b/man/make_shipp_testfiles.Rd new file mode 100644 index 00000000..67c3831e --- /dev/null +++ b/man/make_shipp_testfiles.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/test-helpers.R +\name{make_shipp_testfiles} +\alias{make_shipp_testfiles} +\title{Build JSON from template and a set of params} +\usage{ +make_shipp_testfiles(naomi_output) +} +\arguments{ +\item{naomi_output}{Calibrated naomi output} +} +\value{ +Calibrated naomi output matched to MWI test data on +\code{naomi.resources} to be used to generate the shipp tool. +} +\description{ +Build JSON from template and a set of params +} +\keyword{internal} diff --git a/man/odds.Rd b/man/odds.Rd new file mode 100644 index 00000000..2b7ca872 --- /dev/null +++ b/man/odds.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{odds} +\alias{odds} +\title{Calculate the odds} +\usage{ +odds(p) +} +\arguments{ +\item{p}{Probability in \link{0, 1}} +} +\description{ +Calculate the odds +} diff --git a/man/shipp_adjust_sexbehav_fsw.Rd b/man/shipp_adjust_sexbehav_fsw.Rd new file mode 100644 index 00000000..f98b36b6 --- /dev/null +++ b/man/shipp_adjust_sexbehav_fsw.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{shipp_adjust_sexbehav_fsw} +\alias{shipp_adjust_sexbehav_fsw} +\title{Adjust female sexual behavior risk groups by FSW proportions} +\usage{ +shipp_adjust_sexbehav_fsw(outputs, options, fsw_est) +} +\arguments{ +\item{outputs}{Naomi output.} + +\item{options}{Naomi model options.} + +\item{fsw_est}{5-year estimates of FSW PSEs generated from \code{shipp_disaggregate_fsw()}. + +Estimates are generated for the following groups: +\itemize{ +\item \code{nosex12m}: +\item \code{sexcohab}: +\item \code{sexnonregplus}: +\item \code{sexnonreg}: +\item \code{sexpaid12m}: +\item \code{nosex12m}: +} + +Calculation steps: +\enumerate{ +\item Align admin0/admin1 FSW proportions with SRB SAE estimates. +\item Subtract the proportion of FSW from total high risk female population. +}} +} +\value{ +District level estimates of female sexual risk behaviour groups +} +\description{ +Adjust female sexual behavior risk groups by FSW proportions +} +\keyword{internal} diff --git a/man/shipp_adjust_sexbehav_msm_pwid.Rd b/man/shipp_adjust_sexbehav_msm_pwid.Rd new file mode 100644 index 00000000..467b2ead --- /dev/null +++ b/man/shipp_adjust_sexbehav_msm_pwid.Rd @@ -0,0 +1,57 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{shipp_adjust_sexbehav_msm_pwid} +\alias{shipp_adjust_sexbehav_msm_pwid} +\title{Adjust male sexual behavior risk groups by MSM + PWID proportions} +\usage{ +shipp_adjust_sexbehav_msm_pwid(outputs, options, msm_est, pwid_est) +} +\arguments{ +\item{outputs}{Naomi output.} + +\item{options}{Naomi model options.} + +\item{msm_est}{5-year estimates of MSM PSEs generated from \code{shipp__disaggregate_msm()}.} + +\item{pwid_est}{5-year estimates of MSM PSEs generated from \code{shipp__disaggregate_pwid()}.} +} +\value{ +District level estimates of male sexual risk behaviour groups +} +\description{ +Adjust male sexual behavior risk groups by MSM + PWID proportions +} +\keyword{*} +\keyword{1.} +\keyword{2.} +\keyword{Align} +\keyword{Calculation} +\keyword{Estimates} +\keyword{MSM} +\keyword{PWID} +\keyword{SAE} +\keyword{SRB} +\keyword{Subtracting} +\keyword{`msm`:} +\keyword{`nosex12m`:} +\keyword{`pwid`:} +\keyword{`sexcohab`:} +\keyword{`sexonreg`:} +\keyword{`sexonregplus`:} +\keyword{admin0/admin1} +\keyword{all} +\keyword{and} +\keyword{are} +\keyword{estimates.} +\keyword{following} +\keyword{for} +\keyword{from} +\keyword{generated} +\keyword{groups.} +\keyword{groups:} +\keyword{internal} +\keyword{proportionally} +\keyword{proportions} +\keyword{steps:} +\keyword{the} +\keyword{with} diff --git a/man/shipp_calculate_incidence_female.Rd b/man/shipp_calculate_incidence_female.Rd new file mode 100644 index 00000000..a8352037 --- /dev/null +++ b/man/shipp_calculate_incidence_female.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{shipp_calculate_incidence_female} +\alias{shipp_calculate_incidence_female} +\title{Calculate incidence for female SRB groups.} +\usage{ +shipp_calculate_incidence_female( + naomi_output, + options, + female_srb, + female_logit_prevalence, + survey_year = 2018, + kp_consensus +) +} +\arguments{ +\item{naomi_output}{Naomi output.} + +\item{options}{Naomi model options.} + +\item{female_srb}{FSW adjusted estimates of female sexual risk groups generated by \code{shipp_adjust_sexbehav_fsw()}.} + +\item{female_logit_prevalence}{Risk adjusted estimates of female prevalence in sexual risk groups generated by \code{shipp_calculate_prevalence_female()}.} + +\item{survey_year}{Survey year to sample from the SAE model. Default is 2018. Survey year should be updated to most current household survey in the country - for countries without recent household surveys, leave at 2018 - the spatiotemporal +model of sexual behaviour fitted to all countries has the most data for in roughly 2018.} +} +\value{ +Wide format output required for the SHIPP workbook. +} +\description{ +While maintaining age/sex/district-specific HIV incidence from Naomi, distribute +HIV incidence between our 4 different behavioural groups utilizing IRRs from the +literature +} +\keyword{internal} diff --git a/man/shipp_calculate_incidence_male.Rd b/man/shipp_calculate_incidence_male.Rd new file mode 100644 index 00000000..af301638 --- /dev/null +++ b/man/shipp_calculate_incidence_male.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{shipp_calculate_incidence_male} +\alias{shipp_calculate_incidence_male} +\title{Calculate incidence in high risk male key populations} +\usage{ +shipp_calculate_incidence_male( + naomi_output, + options, + male_srb, + male_logit_prevalence, + survey_year = 2018, + kp_consensus +) +} +\arguments{ +\item{options}{Naomi model options.} + +\item{male_srb}{MSM and PWID adjusted estimated of male SRB groups generated by \code{shipp_adjust_sexbehav_msm_pwid()}.} + +\item{male_logit_prevalence}{Risk adjusted estimates of male HIV prevalence in sexual risk groups generated by \code{shipp_calculate_prevalence_male()}.} + +\item{survey_year}{Survey year to sample from the SAE model. Default is 2018. Survey year should be updated to most current household survey in the country - for countries without recent household surveys, leave at 2018 - the spatiotemporal +model of sexual behaviour fitted to all countries has the most data for in roughly 2018.} + +\item{outputs}{Naomi output.} +} +\value{ +Wide format output required for the SHIPP workbook +} +\description{ +Calculate incidence in high risk male key populations +} diff --git a/man/shipp_calculate_prevalence_female.Rd b/man/shipp_calculate_prevalence_female.Rd new file mode 100644 index 00000000..9ac9e421 --- /dev/null +++ b/man/shipp_calculate_prevalence_female.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{shipp_calculate_prevalence_female} +\alias{shipp_calculate_prevalence_female} +\title{Calculate prevalence for female SRB groups.} +\usage{ +shipp_calculate_prevalence_female( + naomi_output, + options, + fsw_est, + female_srb, + survey_year_sample = 2018 +) +} +\arguments{ +\item{naomi_output}{Naomi output.} + +\item{options}{Naomi model options.} + +\item{fsw_est}{5-year estimates of FSW PSEs generated from \code{shipp_disaggregate_fse()}.} + +\item{female_srb}{FSW adjusted estimates of female SRB groups generated by \code{shipp_adjust_sexbehav_fsw()}} + +\item{survey_year}{Year of survey to sample estimates. + +To calculate district-age-sex-sexual behaviour-specific HIV prevalence, we maintain +HIV prevalence from Naomi for a district-age-sex, but disaggregate to different +risk behaviours using: +(1) HIV prevalence ratios from household surveys for those reporting no +sex vs one cohabiting vs non-regular sexual partner(s), and +(2) a linear regression through admin-1 level estimates of the ratio of KP +prevalence to gen-pop prevalence used to predict an age-district-specific +FSW to general population prevalence ratio.} +} +\value{ +SRB PSEs with logit prevalence estimates. +} +\description{ +Calculate prevalence for female SRB groups. +} +\keyword{internal} diff --git a/man/shipp_calculate_prevalence_male.Rd b/man/shipp_calculate_prevalence_male.Rd new file mode 100644 index 00000000..0ca15b62 --- /dev/null +++ b/man/shipp_calculate_prevalence_male.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{shipp_calculate_prevalence_male} +\alias{shipp_calculate_prevalence_male} +\title{Calculate prevalence for male SRB groups.} +\usage{ +shipp_calculate_prevalence_male( + naomi_output, + areas, + options, + msm_est, + male_srb, + survey_year_sample = 2018 +) +} +\arguments{ +\item{naomi_output}{Naomi output.} + +\item{areas}{Naomi boundary file.} + +\item{options}{Naomi model options.} + +\item{msm_est}{5-year estimates of MSM PSEs generated from \code{shipp__disaggregate_msm()}.} + +\item{male_srb}{MSM and PWID adjusted estimates of male SRB groups generated by \code{shipp_adjust_sexbehav_msm_pwid()}.} + +\item{survey_year}{Year of survey to sample estimates. + +To calculate district-age-sex-sexual behaviour-specific HIV prevalence, we maintain +HIV prevalence from Naomi for a district-age-sex, but disaggregate to different +risk behaviours using: +(1) HIV prevalence ratios from household surveys for those reporting no +sex vs one cohabiting vs non-regular sexual partner(s), and +(2) admin-1 level estimates of the ratio of KP prevalence to gen-pop prevalence +among 15-24 year olds for MSM (due to the young age distribution of MSM) or +among 15-49 year olds for PWID (due to the older age distribution of PWID) +applied to all age groups among MSM and PWID in districts by admin-1 unit.} +} +\value{ +SRB PSEs with logit prevalence estimates. +} +\description{ +Calculate prevalence for male SRB groups. +} +\keyword{internal} diff --git a/man/shipp_disaggregate_fsw.Rd b/man/shipp_disaggregate_fsw.Rd new file mode 100644 index 00000000..3c87611c --- /dev/null +++ b/man/shipp_disaggregate_fsw.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{shipp_disaggregate_fsw} +\alias{shipp_disaggregate_fsw} +\title{Dissagreggate admin1 FSW proportions from Oli's KP model to 5-age groups} +\usage{ +shipp_disaggregate_fsw(outputs, options, naomi_pop, kp_consensus) +} +\arguments{ +\item{outputs}{Naomi output.} + +\item{options}{Naomi model options.} + +\item{naomi_pop}{Naomi population estimates for T2.} + +\item{kp_consensus}{Key pop consensus estimates.} +} +\value{ +District level FSW estimates by 5-year age bands for ages 15-49. +} +\description{ +Dissagreggate admin1 FSW proportions from Oli's KP model to 5-age groups +} +\keyword{internal} diff --git a/man/shipp_disaggregate_msm.Rd b/man/shipp_disaggregate_msm.Rd new file mode 100644 index 00000000..1ba051cb --- /dev/null +++ b/man/shipp_disaggregate_msm.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{shipp_disaggregate_msm} +\alias{shipp_disaggregate_msm} +\title{Disaggregate admin1 MSM proportions from Oli's KP model to 5-age groups} +\usage{ +shipp_disaggregate_msm(outputs, options, naomi_pop, kp_consensus) +} +\arguments{ +\item{outputs}{Naomi output.} + +\item{options}{Naomi model options.} + +\item{naomi_pop}{Naomi population estimates for T2.} + +\item{kp_consensus}{Key pop consensus estimates.} +} +\value{ +District level MSM estimates by 5-year age bands for ages 15-49. +} +\description{ +Disaggregate admin1 MSM proportions from Oli's KP model to 5-age groups +} +\keyword{internal} diff --git a/man/shipp_disaggregate_pwid.Rd b/man/shipp_disaggregate_pwid.Rd new file mode 100644 index 00000000..e40ce820 --- /dev/null +++ b/man/shipp_disaggregate_pwid.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{shipp_disaggregate_pwid} +\alias{shipp_disaggregate_pwid} +\title{Disaggregate admin1 PWID proportions from Oli's KP model to 5-age groups} +\usage{ +shipp_disaggregate_pwid(outputs, options, naomi_pop, kp_consensus) +} +\arguments{ +\item{outputs}{Naomi output.} + +\item{options}{Naomi model options.} + +\item{naomi_pop}{Naomi population estimates for T2.} + +\item{kp_consensus}{Key pop consensus estimates.} +} +\value{ +District level PWID estimates by 5-year age bands for ages 15-49. +} +\description{ +Disaggregate admin1 PWID proportions from Oli's KP model to 5-age groups +} +\keyword{internal} diff --git a/man/shipp_format_naomi.Rd b/man/shipp_format_naomi.Rd new file mode 100644 index 00000000..571629fe --- /dev/null +++ b/man/shipp_format_naomi.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{shipp_format_naomi} +\alias{shipp_format_naomi} +\title{Format naomi outputs for PSE tool} +\usage{ +shipp_format_naomi(outputs, options) +} +\arguments{ +\item{outputs}{Naomi output} + +\item{options}{Naomi model options.} +} +\value{ +Naomi indicators formatted for the SHIPP workbook. +} +\description{ +Format naomi outputs for PSE tool +} +\keyword{internal} diff --git a/man/shipp_generate_risk_populations.Rd b/man/shipp_generate_risk_populations.Rd new file mode 100644 index 00000000..3e7925d2 --- /dev/null +++ b/man/shipp_generate_risk_populations.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/shipp.R +\name{shipp_generate_risk_populations} +\alias{shipp_generate_risk_populations} +\title{Generate outputs to update SHIPP tool.} +\usage{ +shipp_generate_risk_populations(naomi_output, pjnz, survey_year = 2018) +} +\arguments{ +\item{naomi_output}{Path to naomi output (zip file or hintr object).} + +\item{pjnz}{Path to spectrum file.} + +\item{survey_year}{Survey year to sample from the SAE model. Default is 2018. Survey year should be updated to most current household survey in the country - for countries without recent household surveys, leave at 2018 - the spatiotemporal +model of sexual behaviour fitted to all countries has the most data for in roughly 2018.} + +\item{male_srb}{Estimates of male sexual risk groups generated by \code{shipp_adjust_sexbehav_msm_pwid()}} + +\item{male_logit_prevalence}{Risk adjusted estimates of male prevalence in sexual risk groups generated by \code{shipp_calculate_prevalence_male()}} +} +\value{ +Output files to update SHIPP excel workbook. +} +\description{ +Generate outputs to update SHIPP tool. +} +\keyword{internal} diff --git a/man/write_xlsx_sheets.Rd b/man/write_xlsx_sheets.Rd new file mode 100644 index 00000000..f8dafe85 --- /dev/null +++ b/man/write_xlsx_sheets.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{write_xlsx_sheets} +\alias{write_xlsx_sheets} +\title{Write list of data frames into an xlsx file} +\usage{ +write_xlsx_sheets(template, sheets, path) +} +\arguments{ +\item{template}{Path to xlsx file with sheets} + +\item{sheets}{Named list of data frames to write into template. The names +must match the destination sheet in the xlsx} + +\item{path}{Path to output the filled in xlsx} +} +\value{ +Path to complete xlsx file +} +\description{ +This doesn't write colmn headers into the workbook, it expects +that these already exist. +} +\keyword{internal} diff --git a/tests/testthat/helper-testing.R b/tests/testthat/helper-testing.R index 137fda0b..d06ec7fa 100644 --- a/tests/testthat/helper-testing.R +++ b/tests/testthat/helper-testing.R @@ -3,6 +3,8 @@ expect_no_error <- function(object) { } with_mock <- function(..., .parent = parent.frame()) { + ## Don't use this, this should be removed, now testthat have + ## added with_mocked_bindings that should be preferred mockr::with_mock(..., .parent = .parent, .env = "naomi") } diff --git a/tests/testthat/test-downloads.R b/tests/testthat/test-downloads.R index f439f6cf..768f0037 100644 --- a/tests/testthat/test-downloads.R +++ b/tests/testthat/test-downloads.R @@ -192,29 +192,139 @@ test_that("comparison report download can be created", { "Generating comparison report") }) -test_that("AGYW download can be created", { + +test_that("SHIPP download can be created", { + + shipp_output_demo <- make_shipp_testfiles(a_hintr_output_calibrated) + mock_new_simple_progress <- mockery::mock(MockSimpleProgress$new()) - with_mock(new_simple_progress = mock_new_simple_progress, { + + with_mocked_bindings( messages <- naomi_evaluate_promise( - out <- hintr_prepare_agyw_download(a_hintr_output_calibrated, - a_hintr_data$pjnz)) - }) + out <- hintr_prepare_shipp_download(shipp_output_demo, + a_hintr_data$pjnz)), + new_simple_progress = mock_new_simple_progress) + expect_true(file.exists(out$path)) expect_type(out$metadata$description, "character") expect_length(out$metadata$description, 1) expect_equal(out$metadata$areas, "MWI") - read <- readxl::read_xlsx(out$path) - expect_equal(read, - data.frame(x = c(1, 2, 3), y = c(3, 4, 5)), - ignore_attr = TRUE) + outputs_female <- openxlsx::readWorkbook(out$path, sheet = "All outputs - F") + expect_true(nrow(outputs_female) > 10) + outputs_male <- openxlsx::readWorkbook(out$path, sheet = "All outputs - M") + expect_true(nrow(outputs_male) > 10) + naomi_outputs <- openxlsx::readWorkbook(out$path, sheet = "NAOMI outputs") + expect_true(nrow(naomi_outputs) > 4) ## Progress messages printed expect_length(messages$progress, 1) - expect_equal(messages$progress[[1]]$message, "Generating AGYW tool") + expect_equal(messages$progress[[1]]$message, "Generating SHIPP tool") + + # Test shipp workbook with no kp workbook saved into spectrum + risk_prop <- shipp_generate_risk_populations(shipp_output_demo$model_output_path, + a_hintr_data$pjnz) + + expect_equal(risk_prop$meta_consensus, + data.frame(kp = c("FSW", "MSM", "PWID"), + consensus_estimate = NA)) + + # Test shipp workbook with mock workbook saved into spectrum + kp_consensus <- readRDS(file.path("testdata/kp_workbook_spectrum.rds")) + mock_extract_kp_workbook <- mockery::mock(kp_consensus) + mock_new_simple_progress <- mockery::mock(MockSimpleProgress$new()) + + with_mocked_bindings( + risk_prop_scaled <- shipp_generate_risk_populations( + shipp_output_demo$model_output_path, a_hintr_data$pjnz), + new_simple_progress = mock_new_simple_progress, + extract_kp_workbook = mock_extract_kp_workbook + ) + + # Check that consensus estimates extracted and saved out + expect_equal(risk_prop_scaled$meta_consensus, + data.frame(kp = c("FSW", "MSM", "PWID"), + consensus_estimate = c(40000, 35500, 5000))) + + # Test that PSE tool adjusted to KP consensus estimates correctly + model_object <- read_hintr_output(shipp_output_demo$model_output_path) + outputs <- model_object$output_package + options <- outputs$fit$model_options + naomi <- shipp_format_naomi(outputs, options) + + # Naomi population + naomi_pop <- naomi$naomi_long %>% + dplyr::filter(indicator == "population") %>% + dplyr::select(area_id, area_level,sex, age_group, area_level, + spectrum_region_code, population = mean) + + naomi_pop$iso3 <- options$area_scope + + # KP PSEs adjusted to consensus estimates when consensus estimates are + # < 5% of age matched population denominator + fsw_est <- shipp_disaggregate_fsw(outputs, options, naomi_pop, kp_consensus) + pwid_est <- shipp_disaggregate_pwid(outputs, options, naomi_pop, kp_consensus) + msm_est <- shipp_disaggregate_msm(outputs, options, naomi_pop, kp_consensus) + + fsw <- sum(fsw_est$fsw) + pwid <- sum(pwid_est$pwid) + msm <- sum(msm_est$msm) + + # Note that PWID will be 90% of KP workbook consensus estimate due to exclusion + # of female PWID + expect_equal(c(fsw, pwid, msm), c(40000, 4550, 35500)) + + + # KP PSEs **not** adjusted to consensus estimates when consensus estimates are + # > 5% of age matched population denominator + kp_consensus_bad <- readRDS(file.path("testdata/kp_workbook_spectrum_bad.rds")) + mock_extract_kp_workbook <- mockery::mock(kp_consensus_bad) + mock_new_simple_progress <- mockery::mock(MockSimpleProgress$new()) + + with_mocked_bindings( + risk_prop_scaled <- shipp_generate_risk_populations( + shipp_output_demo$model_output_path, a_hintr_data$pjnz), + new_simple_progress = mock_new_simple_progress, + extract_kp_workbook = mock_extract_kp_workbook + ) + + # Check that bad consensus estimates extracted and saved out + expect_equal(risk_prop_scaled$meta_consensus, + data.frame(kp = c("FSW", "MSM", "PWID"), + consensus_estimate = c(260000, 260000, 260000))) + + + # KP PSEs use default proportions from Oli's mode when consensus estimates are + # >= 5% of age matched population denominator + + fsw_est <- shipp_disaggregate_fsw(outputs, options, naomi_pop, kp_consensus_bad) + pwid_est <- shipp_disaggregate_pwid(outputs, options, naomi_pop, kp_consensus_bad) + msm_est <- shipp_disaggregate_msm(outputs, options, naomi_pop, kp_consensus_bad) + + fsw <- sum(fsw_est$fsw) + pwid <- sum(pwid_est$pwid) + msm <- sum(msm_est$msm) + + expect_equal(c(fsw, pwid, msm), c(62306.311, 7398.486, 24794.842)) + +}) + + +test_that("Error thrown when SHIPP resources are out of date", { + + kp_error <- paste0("Available KP PSE estimates for: \n", + "MWI_1_1; MWI_1_2; MWI_1_3", + "\n\n Do not match Naomi estimates for: \n", + "MWI_2_1_demo; MWI_2_2_demo; MWI_2_3_demo; MWI_2_4_demo; MWI_2_5_demo", + "\n\nTo update estimates, please contact Naomi support.") + + expect_error(hintr_prepare_shipp_download(a_hintr_output_calibrated, + a_hintr_data$pjnz), kp_error) + }) + test_that("output description is translated", { text <- build_output_description(a_hintr_options) expect_match( @@ -227,3 +337,10 @@ test_that("output description is translated", { expect_match(text, paste0("Paquet Naomi téléchargée depuis l'application ", "web Naomi\\n\\nPérimètre de zone - MWI\\n.+")) }) + +test_that("failing to write data into xlsx sheet gives a useful error", { + sheets_to_write <- list(x = data.frame(x = c(1, 2, 3))) + dest <- tempfile() + expect_error(write_shipp_workbook(sheets_to_write, dest), + "Failed to build workbook, please contact support: Sheet 'x' does not exist") +}) diff --git a/tests/testthat/testdata/kp_workbook_spectrum.rds b/tests/testthat/testdata/kp_workbook_spectrum.rds new file mode 100644 index 00000000..2477dd72 Binary files /dev/null and b/tests/testthat/testdata/kp_workbook_spectrum.rds differ diff --git a/tests/testthat/testdata/kp_workbook_spectrum_bad.rds b/tests/testthat/testdata/kp_workbook_spectrum_bad.rds new file mode 100644 index 00000000..871fd1f0 Binary files /dev/null and b/tests/testthat/testdata/kp_workbook_spectrum_bad.rds differ diff --git a/vignettes/hiv-prev-workflow.Rmd b/vignettes/hiv-prev-workflow.Rmd new file mode 100644 index 00000000..75c5a84b --- /dev/null +++ b/vignettes/hiv-prev-workflow.Rmd @@ -0,0 +1,189 @@ +--- +title: "SHIPP tool" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{HIV prevention prioritization tool workflow} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, echo = FALSE, warning = FALSE, results = 'asis'} +library(tibble) +library(gt) +library(naomi) +``` + +## Background + +Many HIV prevention programmes aim to reduce new infections but it is infeasible to target all individuals living in locations with moderate to high HIV incidence. In different geographic settings, with different background incidence, sexual behaviours confer different levels of risk of HIV infection. + +Recent HIV programming guidance introduces thresholds for prioritization considering both HIV incidence and sexual behaviours to reach the largest population at risk of HIV. The SHIPP tool provides the “denominator” for HIV prevention categorised by sex, age, geography, and sexual behaviour. + +## Data inputs + +These categorizations are calculated using subnational estimates of HIV prevalence and incidence by age and sex produced by the Naomi model. [Naomi](https://github.com/mrc-ide/naomi) is a small-area estimation model for estimating HIV prevalence and PLHIV, ART coverage, and new HIV infections at a district level by sex and five-year age group. The model combines district-level data about multiple outcomes from several sources in a Bayesian statistical model to produce robust indicators of subnational HIV burden. Naomi is used to update annual estimates of subnational HIV burden as part of the [UNAIDS HIV estimates process](https://www.unaids.org/en/dataanalysis/knowyourresponse/HIVdata_estimates). + +The tool synthesises the most recent estimates of subnational HIV prevalence and incidence with outputs from additional statistical models +describing subnational variation in sexual behaviour and estimates of key populations at an elevated risk of HIV infection: + +* We extend the **[Howes et al.](https://journals.plos.org/globalpublichealth/article?id=10.1371/journal.pgph.0001731)** spatiotemporal multinomial model of sexual behaviour proportions to both men and women aged 15-49. This model utilizes household survey data (Demographic and Health Surveys, Population-based HIV Impact Assessments, Multiple Indicator Cluster Surveys, and other country-specific surveys) to develop spatially smoothed district-level estimates of the proportion of the population in each of three behavioural group (no sex, one regular partner, non-regular partner(s)) described below. The model includes terms for age (5-year age groups), a main spatial effect, a spatial interaction term for 5-year age groups, and survey year (a temporal effect). Men's and women's behavioural data is fitted separately, though all countries are fitted in a single model stratified by sex. The models are fit using the INLA package in R using the multinomial-Poisson transformation. +* **[Stevens et al.](https://www.medrxiv.org/content/10.1101/2022.07.27.22278071v2)** model for population size estimates (PSEs) of key populations in sub-Saharan Africa at an admin-1 level (regional), including female sex workers (FSW), men who have sex with men (MSM), and people who inject drugs (PWID). +* **[Nguyen et al.](https://bmcpublichealth.biomedcentral.com/articles/10.1186/s12889-022-13451-y)**: model for age at sexual debut in sub-Saharan Africa. + +The outputs of these models are stored in an external repository, [naomi-resources](http://github/mrc-ide/naomi.resources) and are updated annually to incorporate newly released survey data and to align with changes in geographic areas in country specific administrative boundaries required for planning. + +In addition to estimates produced by subnational models, the tool incorporates consensus estimates of KP population size and HIV incidence that are developed by national HIV estimates teams as part of the annual UNAIDS HIV estimates process. For more information on this exercise please see [14G Key Population Workbook](https://hivtools.unaids.org/hiv-estimates-training-material-en/). + +## Tool workflow + +This tool is now integrated into the [Naomi web application](https://naomi.unaids.org/login) and is generated after completing the Naomi model as described in [22G Naomi sub-national estimates: Creating subnational HIV estimates](https://hivtools.unaids.org/hiv-estimates-training-material-en/). + +If you are running a Naomi model fit with updated administrative boundaries, you may receive an error that the external database containing the sexual behaviour or KP PSE model is out of date: + +```{r, echo = FALSE, results = 'asis', warning = FALSE} + +cat("Error: Available KP PSE estimates for: \n", + "MWI_1_1; MWI_1_2; MWI_1_3", + "\n\n Do not match Naomi estimates for: \n", + "MWI_1_1xc; MWI_2_25d; MWI_2_3cv", + "\n\nTo update estimates, please contact Naomi support.") +``` + +## SHIPP Tool estimates process + +**1. Estimate key population sizes by district and age**: + +Regional KP proportion estimates from [Stevens et al.](https://www.medrxiv.org/content/10.1101/2022.07.27.22278071v2) are disaggregated by age and district. + + * _Nationally adjusted age distributions for MSM and FSW_: For MSM and FSW, estimates for ages 15-49 are age-disaggregated using South Africa's [Thembisa](https://www.thembisa.org/content/downloadPage/Thembisa4_3) model estimates of the age distribution for FSW and MSM - Gamma(29,9) for FSW and Gamma(25,7) for MSM - then adjusted for the distribution of age at sexual debut by sex and country [(Nguyen et al.)](https://bmcpublichealth.biomedcentral.com/articles/10.1186/s12889-022-13451-y) to account for differences between South Africa and the other countries included in the SHIPP tool. + * _Uniformly adjusted age distribution for PWID_: For age distribution of PWID, we utilize a uniform literature estimate of age distribution from [Hines et al](web) across countries (mean age 29.4, SD 7). + * _Calculating total KP population by age and district_: Country specific regional KP proportions from the _Stevens et al_ model are applied to district level population estimates from Naomi and adjusted by age distribution described above. + * This includes an assumption that a nominal number of PWID (~9%) are female who are removed from the denominator. As such, PWID are assumed as male from this point onwards. + + +**2. Separate general population sexual behaviour groups from KP populations calculated in (1)**: + +```{r, echo = FALSE, results = 'asis', warning = FALSE} + +tibble::tribble( + ~"Category", ~ "HIV related risk", + "No sex", "Not sexually active", + "One regular", "Sexually active, one cohabiting/marital partner", + "Non-regular", "Non-regular sexual partner(s)", + "Key Populations", "FSW (women), MSM and PWID (men)" +) %>% + gt() %>% + tab_header("HIV prevention priority groups") %>% + gt::tab_options( + table.align = "left", + heading.align = "left", + column_labels.font.size = "small", + column_labels.background.color = "grey", + table.font.size = "smaller", + data_row.padding = gt::px(3), + row_group.background.color = "lightgrey" + ) + +``` + +Subtract the proportion of KPs from sexual behaviour groups estimated in _Risher et al_ model: + + * FSWs only subtracted from non-regular partner(s) group of females across ages. + * MSM and PWID subtracted proportionally from all male behaviour categories. + +**3. Estimate HIV prevalence by behaviour** + +HIV prevalence ratios by behaviour group are used to distribute PLHIV between behavioural risk groups. + +* Household survey data is used to estimate HIV prevalence in the no sex, one regular and non-regular partner(s) groups to calculate log odds-ratios for each behavioural category. +* HIV prevalence ratios for KPs are based on the ratio of KP HIV prevalence from the _Stevens et al_ model to HIV prevalence among all women (FSW) or men (MSM and PWID). +* HIV prevalence by behaviour is not explicitly presented – it is used to subtract off population sizes to present the population susceptible to HIV (HIV-negative). + + +_For female KPs, HIV prevalence ratios are derived based on:_ + + * A linear regression through regional (admin-1 level) estimates of the ratio of KP prevalence to general population prevalence used to predict an age-district-specific FSW to general population prevalence ratio. + +_For males KPs, HIV prevalence ratios are derived based on:_ + + * Regional (admin-1 level) estimates of the ratio of KP prevalence to general population prevalence among 15-24 year olds for MSM (due to the young age distribution of MSM) or among 15-49 year olds for PWID (due to the older age distribution of PWID) applied to all age groups among MSM and PWID in districts by region (admin-1 unit). + +**4. Estimate HIV incidence rates and new HIV infections by behaviour** + +While maintaining age/sex/district-specific HIV incidence from Naomi, distribute HIV incidence between our 4 different behavioural groups utilizing incidence rate ratios (IRRs) from the literature: + +* Risk ratios for non-regular sex partner(s) relative to those with a single cohabiting/marital sex partner for females_1,2,3_ and males._1,2,4_ +* For FSW, ratio of HIV incidence among women in key populations vs general population women derived based on HIV incidence category in district and tiered risk ratio._5_ +* For MSM and PWID, using regional (admin-1) KP prevalence estimates relative to general population prevalence_6_ and estimates from systematic review & meta-regression._7_ +* National KP new infections are scaled to consensus estimates of KP HIV incidence maintaining the relative district-level proportions associated with the above two bullet points. As such, the precise risk ratios for KPs will not necessarily match the risk ratios listed in the table below. + + +Number of new infections by sexual behaviour group is derived by multiplying these estimated HIV incidence rates by behaviour times the population sizes of HIV-negative individuals by behaviour, in each of the 5-year age groups + +```{r, echo = FALSE, results = 'asis', warning = FALSE} + +tibble::tribble( + ~"Category", ~ "Females", ~"Males" , + "No sex", "0", "0", + "One regular", "1 (reference category)", "1 (reference category)", + "Non-regular", + "_Aged 15-24_: 1.72
Aged 25-49: 2.1 _1,2,3_", + "_Aged 15-24_: 1.89
Aged 25-49: 2.1 _1,2,4_", + "Key Populations", + "**FSW**:
Very high: 3
High: 4
Moderate: 7
Low: 11
Very low: 17 _5_", + "**MSM**: 2.5-250_6,7_
**PWID** 2.5-55 _6_" +) %>% + gt() %>% + gt::tab_options( + table.align = "left", + heading.align = "left", + column_labels.font.size = "small", + column_labels.background.color = "grey", + table.font.size = "smaller", + data_row.padding = gt::px(3), + row_group.background.color = "lightgrey" + ) %>% + fmt_markdown(columns = everything()) + +``` + + + + + + 1 ALPHA Network pooled analysis (Slaymaker et al. CROI 2020) + + 2 [Jia et al. 2022](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8743366/) + + 3 [Ssempijja et al. 2022](https://journals.lww.com/jaids/fulltext/2022/07010/high_rates_of_pre_exposure_prophylaxis_eligibility.7.aspx). + + 4 [Hoffman et al. 2022](https://journals.lww.com/jaids/Fulltext/2022/06001/Implementing_PrEP_Services_in_Diverse_Health_Care.15.aspx) + + 5 [Jones et al. 2023](https://www.medrxiv.org/content/10.1101/2023.10.17.23297108v2) + + 6 [Stevens et al. 2023](https://www.medrxiv.org/content/10.1101/2022.07.27.22278071v2) + + 7 [Stannah et al. 2022](https://www.medrxiv.org/content/10.1101/2022.11.14.22282329v1) + + + +## Limitations of tool + +Risk behaviour population size estimates at a district level have a high degree of uncertainty, which is not captured in the current version of the tool. + +There is uncertainty in: + +* KP sizes and their age and geographical disaggregation +* Small area estimates based on survey data +* Behavioural reports in surveys +* Incidence rate ratios by behaviour +* Naomi estimates + +As such, SHIPP tool estimates should be considered indicative rather than exact. + + + +## Usage of tool outputs + +Recommendations for usage of tool outputs for prioritising groups for HIV prevention can be found [on the UNAIDS website](https://hivtools.unaids.org/pse/). +