Skip to content

Commit

Permalink
Merge pull request #43 from mrc-ide/spectrum-art-alloc
Browse files Browse the repository at this point in the history 
implement Spectrum ART allocation
  • Loading branch information
@jeffeaton
jeffeaton authored Dec 6, 2024
2 parents 1f5b8d7 + aa55658 commit af6f0a4
Show file tree
Hide file tree
Showing 6 changed files with 156 additions and 27 deletions.
1 change: 1 addition & 0 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
Package: eppasm
Title: Age-structured EPP Model for HIV Epidemic Estimates
Version: 0.8.1
Version: 0.7.7
Authors@R: person("Jeff", "Eaton", email = "[email protected]", role = c("aut", "cre"))
Description: What the package does (one paragraph).
Expand Down
25 changes: 25 additions & 0 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,28 @@
## 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.

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.

* 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.7

* Update to use full names for R internal functions e.g. `Rf_allocVector` instead of `allocVector`. Shorthand names are no longer allowed in R v4.5.0. See Nov 10th news https://developer.r-project.org/blosxom.cgi/R-devel
Expand Down
48 changes: 33 additions & 15 deletions R/eppasm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down
25 changes: 25 additions & 0 deletions R/read-spectrum-files.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -576,6 +576,29 @@ read_hivproj_param <- function(pjnz, use_ep5=FALSE){
artelig_specpop <- stats::setNames(dpsub("<PopsEligTreat MV>", 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 ("<AdultARTAdjFactor>")
## * Enabled / disabled by checkbox flag ("<AdultARTAdjFactorFlag>")
## * 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("<AdultARTAdjFactorFlag>") &&
dpsub("<AdultARTAdjFactorFlag>", 2, 4) == 1) {

adult_artadj_factor <- sapply(dpsub("<AdultARTAdjFactor>", 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("<NewARTPatAllocationMethod MV2>"))
art_alloc_method <- as.integer(dpsub("<NewARTPatAllocationMethod MV2>", 2, 4))
else
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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,
Expand Down
6 changes: 4 additions & 2 deletions R/spectrum.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)),
Expand Down
78 changes: 68 additions & 10 deletions src/eppasm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down Expand Up @@ -661,7 +662,7 @@ extern "C" {
}
}
}

// add new infections to HIV population
for(int g = 0; g < NG; g++){
int a = 0;
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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])
Expand Down

0 comments on commit af6f0a4

Please sign in to comment.