From 8a1ce3638f74b2d8cfc325f4fa51ec040157ae44 Mon Sep 17 00:00:00 2001 From: jeffeaton Date: Sat, 20 Jul 2024 21:17:13 -0400 Subject: [PATCH 1/3] implement Spectrum ART allocation --- DESCRIPTION | 2 +- NEWS.md | 16 +++++++++++ R/eppasm.R | 48 +++++++++++++++++++++---------- src/eppasm.cpp | 78 +++++++++++++++++++++++++++++++++++++++++++------- 4 files changed, 118 insertions(+), 26 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index b093a0a..0c8d3c6 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: eppasm Title: Age-structured EPP Model for HIV Epidemic Estimates -Version: 0.7.6 +Version: 0.8.0 Authors@R: person("Jeff", "Eaton", email = "jeffrey.eaton@imperial.ac.uk", role = c("aut", "cre")) Description: What the package does (one paragraph). Depends: R (>= 3.1.0), diff --git a/NEWS.md b/NEWS.md index ec1b8ee..7c91cd0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,19 @@ +## eppasm 0.8.0 + +* Implement Spectrum ART allocation. + + There has been a longstanding discrepancy betweeen EPP-ASM and Spectrum in ART allocation. + For ART allocation by 'expected mortality', EPP-ASM allocated according to mortality by CD4 + and age. + + Spectrum allocates ART in a two step process: first, ART is allocated by CD4 category based + on the 'expected mortality' and 'proportional to eligibility' weight. Second, within + CD4 categories, ART is allocated by age solely proportional to number in each age + group (propotional to eligibility). + + This has modest overall difference, but was a source of numerical differences between + Spectrum and EPP-ASM. + ## eppasm 0.7.6 * Update internal data country ISO3 list to contain St. Kitts & Nevis and Dominica diff --git a/R/eppasm.R b/R/eppasm.R index f9942d6..13f3aae 100644 --- a/R/eppasm.R +++ b/R/eppasm.R @@ -329,23 +329,41 @@ simmod.specfp <- function(fp, VERSION="C", ...) { } } else { + + ## ## Old EPP-ASM implementation + ## + ## Applied 'expected mortality' weight within each cd4-age mortality strata. + ## Different from Spectrum -- see Spectrum implementation below + ## + ## expect.mort.weight <- sweep(fp$cd4_mort[, h.age15plus.idx,], 3, + ## colSums(art15plus.elig * fp$cd4_mort[, h.age15plus.idx,],,2), "/") + ## artinit.weight <- sweep(fp$art_alloc_mxweight * expect.mort.weight, 3, (1 - fp$art_alloc_mxweight)/colSums(art15plus.elig,,2), "+") + ## artinit <- pmin(sweep(artinit.weight * art15plus.elig, 3, art15plus.inits, "*"), + ## art15plus.elig) + + ## Spectrum ART initiation is 2-step process + ## 1. Allocate by CD4 category (weighted by 'eligible' and 'expected mortality') + ## 2. Allocate by age groups (weighted only by eligibility) + + ## First step: allocate initiation by CD4 category (_hm) + artelig_hm <- apply(art15plus.elig, c(1, 3), sum) + expected_deaths_hm <- apply(art15plus.elig * fp$cd4_mort[, h.age15plus.idx,], c(1, 3), sum) + expected.mort.weight_hm <- sweep(expected_deaths_hm, 2, colSums(expected_deaths_hm), "/") + artelig.weight_hm <- sweep(artelig_hm, 2, colSums(artelig_hm), "/") - expect.mort.weight <- sweep(fp$cd4_mort[, h.age15plus.idx,], 3, - colSums(art15plus.elig * fp$cd4_mort[, h.age15plus.idx,],,2), "/") - artinit.weight <- sweep(fp$art_alloc_mxweight * expect.mort.weight, 3, (1 - fp$art_alloc_mxweight)/colSums(art15plus.elig,,2), "+") - artinit <- pmin(sweep(artinit.weight * art15plus.elig, 3, art15plus.inits, "*"), - art15plus.elig) - - ## Allocation by average mortality across CD4, trying to match Spectrum - ## artelig_by_cd4 <- apply(art15plus.elig, c(1, 3), sum) - ## expectmort_by_cd4 <- apply(art15plus.elig * fp$cd4_mort[, h.age15plus.idx,], c(1, 3), sum) - - ## artinit_dist <- fp$art_alloc_mxweight * sweep(artelig_by_cd4, 2, colSums(artelig_by_cd4), "/") + - ## (1 - fp$art_alloc_mxweight) * sweep(expectmort_by_cd4, 2, colSums(expectmort_by_cd4), "/") + artinit.weight_hm <- fp$art_alloc_mxweight * expected.mort.weight_hm + + (1.0 - fp$art_alloc_mxweight) * artelig.weight_hm + + artinit_hm <- sweep(artinit.weight_hm, 2, art15plus.inits, "*") + + ## Second step: within each CD4 category, allocate initiation + ## proportionally by age - ## artinit_prob <- sweep(artinit_dist, 2, art15plus.inits, "*") / artelig_by_cd4 - ## artinit <- sweep(art15plus.elig, c(1, 3), artinit_prob, "*") - ## artinit <- pmin(artinit, art15plus.elig, na.rm=TRUE) + ## Proportion initiating in each sex x CD4 category + artinit_prob <- artinit_hm / artelig_hm + artinit <- sweep(art15plus.elig, c(1, 3), artinit_prob, "*") + + artinit <- pmin(artinit, art15plus.elig, na.rm=TRUE) } } else { diff --git a/src/eppasm.cpp b/src/eppasm.cpp index 55906f3..2c18bf3 100644 --- a/src/eppasm.cpp +++ b/src/eppasm.cpp @@ -606,7 +606,8 @@ extern "C" { double artpop_hahm = 0.0; for(int hu = 0; hu < hTS; hu++) artpop_hahm += artpop[t][g][ha][hm][hu]; - cd4mx_scale = hivpop[t][g][ha][hm] / (hivpop[t][g][ha][hm] + artpop_hahm); + cd4mx_scale = (hivpop[t][g][ha][hm] + artpop_hahm) > 0 ? + hivpop[t][g][ha][hm] / (hivpop[t][g][ha][hm] + artpop_hahm) : 1.0; } double deaths = cd4mx_scale * cd4_mort[g][ha][hm] * hivpop[t][g][ha][hm]; @@ -661,7 +662,7 @@ extern "C" { } } } - + // add new infections to HIV population for(int g = 0; g < NG; g++){ int a = 0; @@ -726,20 +727,57 @@ extern "C" { // ART initiation for(int g = 0; g < NG; g++){ - double artelig_hahm[hAG_15PLUS][hDS], Xart_15plus = 0.0, Xartelig_15plus = 0.0, expect_mort_artelig15plus = 0.0; + + // Spectrum ART allocation is 2-step process + // 1. Allocate by CD4 category (weighted by 'eligible' and 'expected mortality') + // 2. Allocate by age groups (weighted only by eligibility) + // + // The first step: allocate initiation by CD4 category (_hm) requires + // tabulating number eligible and expected mortality within each CD4 + // category (aggregated over all ages). + + double Xart_15plus = 0.0; // Total currently on ART + + double artelig_hahm[hAG_15PLUS][hDS]; + double artelig_hm[hDS]; + double Xartelig_15plus = 0.0; + + double expect_mort_artelig_hm[hDS]; + double expect_mort_artelig15plus = 0.0; + + // Initialise to zero + memset(artelig_hm, 0, hDS * sizeof(double)); + memset(expect_mort_artelig_hm, 0, hDS * sizeof(double)); + for(int ha = hIDX_15PLUS; ha < hAG; ha++){ for(int hm = everARTelig_idx; hm < hDS; hm++){ + if(hm >= anyelig_idx){ - double prop_elig = (hm >= cd4elig_idx) ? 1.0 : (hm >= hIDX_CD4_350) ? 1.0 - (1.0-specpop_percelig[t])*(1.0-who34percelig) : specpop_percelig[t]; - Xartelig_15plus += artelig_hahm[ha-hIDX_15PLUS][hm] = prop_elig * hivpop[t][g][ha][hm] ; - expect_mort_artelig15plus += cd4_mort[g][ha][hm] * artelig_hahm[ha-hIDX_15PLUS][hm]; + + // Specify proportion eligibly + double prop_elig = (hm >= cd4elig_idx) ? 1.0 : + (hm >= hIDX_CD4_350) ? + 1.0 - (1.0-specpop_percelig[t]) * (1.0-who34percelig) : + specpop_percelig[t]; + + double artelig_hahm_tmp = prop_elig * hivpop[t][g][ha][hm]; + artelig_hahm[ha-hIDX_15PLUS][hm] = artelig_hahm_tmp; + artelig_hm[hm] += artelig_hahm_tmp; + Xartelig_15plus += artelig_hahm_tmp; + + double expect_mort_hahm = cd4_mort[g][ha][hm] * artelig_hahm_tmp; + expect_mort_artelig_hm[hm] += expect_mort_hahm; + expect_mort_artelig15plus += expect_mort_hahm; } - for(int hu = 0; hu < hTS; hu++) + + for(int hu = 0; hu < hTS; hu++) { Xart_15plus += artpop[t][g][ha][hm][hu] + DT * gradART[g][ha][hm][hu]; + } } // if pw_artelig, add pregnant women to artelig_hahm population if(g == FEMALE && pw_artelig[t] > 0 && ha < hAG_FERT){ + double frr_pop_ha = 0; for(int a = hAG_START[ha]; a < hAG_START[ha]+hAG_SPAN[ha]; a++) frr_pop_ha += pop[t][HIVN][g][a]; // add HIV- population @@ -748,11 +786,17 @@ extern "C" { for(int hu = 0; hu < hTS; hu++) frr_pop_ha += frr_art[t][ha-hIDX_FERT][hm][hu] * artpop[t][g][ha][hm][hu]; } + + // Add pregnant women in CD4 categories before all eligible for(int hm = anyelig_idx; hm < cd4elig_idx; hm++){ double pw_elig_hahm = births_by_ha[ha-hIDX_FERT] * frr_cd4[t][ha-hIDX_FERT][hm] * hivpop[t][g][ha][hm] / frr_pop_ha; artelig_hahm[ha-hIDX_15PLUS][hm] += pw_elig_hahm; + artelig_hm[hm] += pw_elig_hahm; Xartelig_15plus += pw_elig_hahm; - expect_mort_artelig15plus += cd4_mort[g][ha][hm] * pw_elig_hahm; + + double pw_expect_mort_hahm = cd4_mort[g][ha][hm] * pw_elig_hahm; + expect_mort_artelig_hm[hm] += pw_expect_mort_hahm; + expect_mort_artelig15plus += pw_expect_mort_hahm; } } } // loop over ha @@ -862,10 +906,24 @@ extern "C" { } else { // Use mixture of eligibility and expected mortality for initiation distribution + // ART allocation step 1: allocate by CD4 stage + double artinit_hm[hDS]; + for(int hm = anyelig_idx; hm < hDS; hm++){ + artinit_hm[hm] = artinit_hts * + ( (1.0 - art_alloc_mxweight) * artelig_hm[hm] / Xartelig_15plus + + art_alloc_mxweight * expect_mort_artelig_hm[hm] / expect_mort_artelig15plus); + } + for(int ha = hIDX_15PLUS; ha < hAG; ha++) for(int hm = anyelig_idx; hm < hDS; hm++){ - if (Xartelig_15plus > 0.0) { - double artinit_hahm = artinit_hts * artelig_hahm[ha-hIDX_15PLUS][hm] * ((1.0 - art_alloc_mxweight)/Xartelig_15plus + art_alloc_mxweight * cd4_mort[g][ha][hm] / expect_mort_artelig15plus); + + // ART allocation step 2: within CD4 category, allocate + // by age proportional to eligibility + if (artelig_hm[hm] > 0.0) { + + double artinit_hahm = artinit_hm[hm] * + artelig_hahm[ha-hIDX_15PLUS][hm] / artelig_hm[hm]; + if(artinit_hahm > artelig_hahm[ha-hIDX_15PLUS][hm]) artinit_hahm = artelig_hahm[ha-hIDX_15PLUS][hm]; if(artinit_hahm > hivpop[t][g][ha][hm] + DT * grad[g][ha][hm]) From 538c864c7a75c67c783905e9d8f20a78214926d5 Mon Sep 17 00:00:00 2001 From: jeffeaton Date: Thu, 25 Jul 2024 17:39:33 +0100 Subject: [PATCH 2/3] adjust art_dropout to annualised rate --- NEWS.md | 4 ++++ R/spectrum.R | 6 ++++-- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index 7c91cd0..d623d9a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -13,6 +13,10 @@ This has modest overall difference, but was a source of numerical differences between Spectrum and EPP-ASM. + +* Patch ART dropout implementation. Spectrum converts input ART dropout percent to an + annual rate using [dropout rate] = -log(1.0 - [input percent]). + ## eppasm 0.7.6 diff --git a/R/spectrum.R b/R/spectrum.R index c37e3f6..9d1a0b7 100644 --- a/R/spectrum.R +++ b/R/spectrum.R @@ -220,8 +220,10 @@ create_spectrum_fixpar <- function(projp, demp, hiv_steps_per_year = 10L, proj_s } else { fp$art_dropout_recover_cd4 <- art_dropout_recover_cd4 } - - fp$art_dropout <- projp$art_dropout[as.character(proj_start:proj_end)]/100 + + ## Convert input percent dropout in 12 months to an annual rate (Rob Glaubius email 25 July 2024) + fp$art_dropout <- -log(1.0 - projp$art_dropout[as.character(proj_start:proj_end)]/100) + fp$median_cd4init <- projp$median_cd4init[as.character(proj_start:proj_end)] fp$med_cd4init_input <- as.integer(fp$median_cd4init > 0) fp$med_cd4init_cat <- replace(findInterval(-fp$median_cd4init, - c(1000, 500, 350, 250, 200, 100, 50)), From 3e01d9210d0db2bb66fd9505477875ada777f7a0 Mon Sep 17 00:00:00 2001 From: jeffeaton Date: Wed, 16 Oct 2024 23:01:56 +0100 Subject: [PATCH 3/3] implement Spectrum Adult ART adjustment --- DESCRIPTION | 2 +- NEWS.md | 6 ++++++ R/read-spectrum-files.R | 25 +++++++++++++++++++++++++ 3 files changed, 32 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0c8d3c6..cc95c66 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: eppasm Title: Age-structured EPP Model for HIV Epidemic Estimates -Version: 0.8.0 +Version: 0.8.1 Authors@R: person("Jeff", "Eaton", email = "jeffrey.eaton@imperial.ac.uk", role = c("aut", "cre")) Description: What the package does (one paragraph). Depends: R (>= 3.1.0), diff --git a/NEWS.md b/NEWS.md index d623d9a..b14a1cb 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +## epapsm 0.8.1 + +* Implement Spectrum Adult ART scalar adjustment. This is a user input that + allows the input number on ART to be adjusted by a scalar to account for + over/under-reporting of treatment numbers. + ## eppasm 0.8.0 * Implement Spectrum ART allocation. diff --git a/R/read-spectrum-files.R b/R/read-spectrum-files.R index e6677c4..a90f3da 100644 --- a/R/read-spectrum-files.R +++ b/R/read-spectrum-files.R @@ -576,6 +576,29 @@ read_hivproj_param <- function(pjnz, use_ep5=FALSE){ artelig_specpop <- stats::setNames(dpsub("", 3:9, 2:6), c("description", "pop", "elig", "percent", "year")) } + ## # 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 ("") + ## * Scaling factor only applies to number inputs, not percentages (John Stover email, 20 Feb 2023) + ## -> Even if scaling factor specified in a year with percentage input, ignore it. + ## + + if (exists_dptag("") && + dpsub("", 2, 4) == 1) { + + adult_artadj_factor <- sapply(dpsub("", 3:4, timedat.idx), as.numeric) + + ## Only apply if is number (! is percentage) + adult_artadj_factor <- adult_artadj_factor ^ as.numeric(!art15plus_numperc) + + art15plus_num <- art15plus_num * adult_artadj_factor + } else { + adult_artadj_factor <- array(1.0, dim(art15plus_num)) + } + + if(exists_dptag("")) art_alloc_method <- as.integer(dpsub("", 2, 4)) else @@ -610,6 +633,7 @@ read_hivproj_param <- function(pjnz, use_ep5=FALSE){ dimnames(art15plus_numperc) <- list(sex=c("Male", "Female"), year=proj.years) dimnames(art15plus_num) <- list(sex=c("Male", "Female"), year=proj.years) + dimnames(adult_artadj_factor) <- list(sex=c("Male", "Female"), year=proj.years) artelig_specpop$pop <- c("PW", "TBHIV", "DC", "FSW", "MSM", "IDU", "OTHER") artelig_specpop$elig <- as.logical(as.integer(artelig_specpop$elig)) @@ -738,6 +762,7 @@ read_hivproj_param <- function(pjnz, use_ep5=FALSE){ "artmx_timerr" = artmx_timerr, "art15plus_numperc" = art15plus_numperc, "art15plus_num" = art15plus_num, + "adult_artadj_factor" = adult_artadj_factor, "art15plus_eligthresh" = art15plus_eligthresh, "artelig_specpop" = artelig_specpop, "median_cd4init" = median_cd4init,