Merge pull request #45 from mrc-ide/excess-nonaids-mx

Excess nonaids mx
mrc-ide · Dec 6, 2024 · 7c468a7 · 7c468a7
2 parents af6f0a4 + c17fe68
commit 7c468a7
Show file tree

Hide file tree

Showing 6 changed files with 270 additions and 30 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,6 @@
 Package: eppasm
 Title: Age-structured EPP Model for HIV Epidemic Estimates
-Version: 0.8.1
-Version: 0.7.7
+Version: 0.8.2
 Authors@R: person("Jeff", "Eaton", email = "[email protected]", role = c("aut", "cre"))
 Description: What the package does (one paragraph).
 Depends: R (>= 3.1.0),

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,41 @@
+## eppasm 0.8.2
+
+Implement new excess non-AIDS mortality among PLHIV implemented in 
+Spectrum 6.38 for 2025 UNAIDS HIV estimates.
+
+Spectrum 6.38 implements rates of non-AIDS excess mortality by sex, 
+age group, CD4 category and ART status. By default these mortality rates 
+are applied in concentrated epidemic countries and defaulted to zero in
+African HIV epidemic settings.
+
+Spectrum outputs for non-AIDS deaths among PLHIV in AIM include non-AIDS 
+excess deaths. EPP-ASM deaths outputs are updated accordingly.
+
+To ensure backward compatibility, excess non-AIDS mortality are initialized 
+to 0.0 and replaced with values read from Spectrum if the relevant values
+exist in the .DP file.
+
+Specific changes:
+* `read_hivproj_param()` is updated to read input values for excess non-AIDS
+  mortality from Spectrum .DP file, stored in tag `"<AdultNonAIDSExcessMort MV>"`.
+
+* `read_hivproj_output()` updated to read additiounal deaths outputs:
+   - Non-AIDS excess deaths among PLHIV from the tag `"<NonAIDSExcessDeathsSingleAge MV>"`,
+     and output in the array `specres$nonaids_excess_deaths`.
+   - AIDS deaths by ART status are now output by single-year age in outputs `aidsdeaths_art1`
+     and `aidsdeaths_noart1`. (Previously only output by 5-year age group.
+   - **BREAKING CHANGE:** Previous outputs for AIDS deaths by ART status in 5-year age groups
+     `aidsdeaths_art` and `aidsdeaths_noart` have added a suffix `5` to avoid ambiguity:
+	 `aidsdeaths_art5` and `aidsdeaths_noart5`. Previous code that used `aidsdeaths_art` and
+	 `aidsdeaths_noart` will need to be updated.
+
+* EPP-ASM code updated to incorporate excess non-AIDS mortality among PLHIV. 
+  Internally new model parameters arrays `cd4_nonaids_excess_mort` and 
+  `art_nonaids_excess_mort` follow the same dimensions and stratification 
+  of `cd4_mort` and `art_mort` arrays. These represent expansions of the 
+  arrays saved in Spectrum-AIM, to all age groups and ART duration categories, 
+  consistent with handling of the `cd4_mort` and `art_mort` arrays.
+
 ## epapsm 0.8.1
 
 * Implement Spectrum Adult ART scalar adjustment. This is a user input that 

diff --git a/R/eppasm.R b/R/eppasm.R
@@ -60,6 +60,16 @@ simmod.specfp <- function(fp, VERSION="C", ...) {
   hivdeaths <- array(0, c(pAG, NG, PROJ_YEARS))
   natdeaths <- array(0, c(pAG, NG, PROJ_YEARS))
 
+  excessnonaidsdeaths <- array(0.0, c(pAG, NG, PROJ_YEARS))
+
+  aidsdeaths_noart <- array(0.0, c(hDS, hAG, NG, PROJ_YEARS))
+  natdeaths_noart <- array(0.0, c(hDS, hAG, NG, PROJ_YEARS))
+  excessnonaidsdeaths_noart <- array(0.0, c(hDS, hAG, NG, PROJ_YEARS))
+
+  aidsdeaths_art <- array(0.0, c(hTS, hDS, hAG, NG, PROJ_YEARS))
+  natdeaths_art <- array(0.0, c(hTS, hDS, hAG, NG, PROJ_YEARS))
+  excessnonaidsdeaths_art <- array(0.0, c(hTS, hDS, hAG, NG, PROJ_YEARS))
+
   popadj.prob <- array(0, c(pAG, NG, PROJ_YEARS))
 
   if(fp$eppmod != "directincid_ann") {
@@ -126,9 +136,12 @@ simmod.specfp <- function(fp, VERSION="C", ...) {
     pop[,,,i] <- pop[,,,i] - deaths
     natdeaths[,,i] <- rowSums(deaths,,2)
 
+    natdeaths_noart[,,,i] <- sweep(hivpop[,,,i], 2:3, 1.0 - hiv.sx.prob, "*")
     hivpop[,,,i] <- sweep(hivpop[,,,i], 2:3, hiv.sx.prob, "*")
-    if(i > fp$tARTstart)
+    if(i > fp$tARTstart) {
+      natdeaths_art[,,,,i] <- sweep(artpop[,,,,i], 3:4, 1.0 - hiv.sx.prob, "*")
       artpop[,,,,i] <- sweep(artpop[,,,,i], 3:4, hiv.sx.prob, "*")
+    }
 
     if (fp$projection_period == "midyear") {
 
@@ -217,18 +230,20 @@ simmod.specfp <- function(fp, VERSION="C", ...) {
         cd4mx_scale <- hivpop[,,,i] / (hivpop[,,,i] + colSums(artpop[,,,,i]))
         cd4mx_scale[!is.finite(cd4mx_scale)] <- 1.0
         cd4_mort_ts <- fp$cd4_mort * cd4mx_scale
-      } else
+      } else {
         cd4_mort_ts <- fp$cd4_mort
+      }
 
       grad <- grad - cd4_mort_ts * hivpop[,,,i]              # HIV mortality, untreated
+      hivdeaths.ts <- DT * colSums(cd4_mort_ts * hivpop[,,,i])
+      aidsdeaths_noart[,,,i] <- aidsdeaths_noart[,,,i] + DT * cd4_mort_ts * hivpop[,,,i]
+
+      ## Non-AIDS excess mortality
+      grad <- grad - fp$cd4_nonaids_excess_mort * hivpop[,,,i]
+      nonaids_excess.ts <- DT * colSums(fp$cd4_nonaids_excess_mort * hivpop[,,,i])
+      excessnonaidsdeaths_noart[,,,i] <- excessnonaidsdeaths_noart[,,,i] + DT * fp$cd4_nonaids_excess_mort * hivpop[,,,i]
 
-      ## Remove hivdeaths from pop
-      hivdeaths.ts <- DT*(colSums(cd4_mort_ts * hivpop[,,,i]) + colSums(fp$art_mort * fp$artmx_timerr[ , i] * artpop[,,,,i],,2))
-      calc.agdist <- function(x) {d <- x/rep(ctapply(x, ag.idx, sum), h.ag.span); d[is.na(d)] <- 0; d}
-      hivdeaths_p.ts <- apply(hivdeaths.ts, 2, rep, h.ag.span) * apply(pop[,,hivp.idx,i], 2, calc.agdist)  # HIV deaths by single-year age
-      pop[,,2,i] <- pop[,,2,i] - hivdeaths_p.ts
-      hivdeaths[,,i] <- hivdeaths[,,i] + hivdeaths_p.ts
-
+
       ## ART initiation
       if(i >= fp$tARTstart) {
 
@@ -241,6 +256,18 @@ simmod.specfp <- function(fp, VERSION="C", ...) {
 
         gradART <- gradART - fp$art_mort * fp$artmx_timerr[ , i] * artpop[,,,,i]   # ART mortality
 
+        hivdeaths.ts <- hivdeaths.ts +
+          DT * colSums(fp$art_mort * fp$artmx_timerr[ , i] * artpop[,,,,i],,2)
+        aidsdeaths_art[,,,,i] <- aidsdeaths_art[,,,,i] +
+          DT * fp$art_mort * fp$artmx_timerr[ , i] * artpop[,,,,i]
+
+        ## Non-AIDS excess mortality
+        gradART <- gradART - fp$art_nonaids_excess_mort * artpop[,,,,i]
+        nonaids_excess.ts <- nonaids_excess.ts +
+          DT * colSums(fp$art_nonaids_excess_mort * artpop[,,,,i],,2)
+        excessnonaidsdeaths_art[,,,,i] <- excessnonaidsdeaths_art[,,,,i] +
+          DT * fp$art_nonaids_excess_mort * artpop[,,,,i]
+
 
         ## ART dropout
         ## remove proportion from all adult ART groups back to untreated pop
@@ -404,7 +431,18 @@ simmod.specfp <- function(fp, VERSION="C", ...) {
       }
 
       hivpop[,,,i] <- hivpop[,,,i] + DT * grad
-    }    
+
+
+      ## Remove hivdeaths from pop
+      calc.agdist <- function(x) {d <- x/rep(ctapply(x, ag.idx, sum), h.ag.span); d[is.na(d)] <- 0; d}
+      hivdeaths_p.ts <- apply(hivdeaths.ts, 2, rep, h.ag.span) * apply(pop[,,hivp.idx,i], 2, calc.agdist)  # HIV deaths by single-year age
+      pop[,,2,i] <- pop[,,2,i] - hivdeaths_p.ts
+      hivdeaths[,,i] <- hivdeaths[,,i] + hivdeaths_p.ts
+
+      nonaids_excess_p.ts <- apply(nonaids_excess.ts, 2, rep, h.ag.span) * apply(pop[,,hivp.idx,i], 2, calc.agdist)  # Non-AIDS excess deaths by single-year age
+      pop[,,2,i] <- pop[,,2,i] - nonaids_excess_p.ts
+      excessnonaidsdeaths[,,i] <- excessnonaidsdeaths[,,i] + nonaids_excess_p.ts
+    }
 
     ## ## Code for calculating new infections once per year to match prevalence (like Spectrum)
     ## ## incidence
@@ -432,7 +470,7 @@ simmod.specfp <- function(fp, VERSION="C", ...) {
       pop[,,,i] <- sweep(pop[,,,i], 1:2, mr.prob, "*")
 
       hivpop[,,,i] <- sweep(hivpop[,,,i], 2:3, hiv.mr.prob, "*")
-      if(i > fp$tARTstart)
+      if(i >= fp$tARTstart)
         artpop[,,,,i] <- sweep(artpop[,,,,i], 3:4, hiv.mr.prob, "*")
     }
 
@@ -481,6 +519,16 @@ simmod.specfp <- function(fp, VERSION="C", ...) {
   attr(pop, "hivdeaths") <- hivdeaths
   attr(pop, "natdeaths") <- natdeaths
 
+  attr(pop, "excessnonaidsdeaths") <- excessnonaidsdeaths
+
+  attr(pop, "aidsdeaths_noart") <- aidsdeaths_noart
+  attr(pop, "natdeaths_noart") <- natdeaths_noart
+  attr(pop, "excessnonaidsdeaths_noart") <- excessnonaidsdeaths_noart
+
+  attr(pop, "aidsdeaths_art") <- aidsdeaths_art
+  attr(pop, "natdeaths_art") <- natdeaths_art
+  attr(pop, "excessnonaidsdeaths_art") <- excessnonaidsdeaths_art
+
   attr(pop, "popadjust") <- popadj.prob
 
   attr(pop, "pregprevlag") <- pregprevlag

diff --git a/R/read-spectrum-files.R b/R/read-spectrum-files.R
@@ -238,10 +238,10 @@ read_hivproj_output <- function(pjnz, single.age=TRUE){
 
   aidsdeaths5 <- array(sapply(aidsdeaths5, as.numeric), lengths(dn5), dn5)
 
-  aidsdeaths_art <- array(sapply(dpsub("<AIDSDeathsART MV2>", c(4:20, 22:38), timedat.idx), as.numeric),
-                          lengths(dn5), dn5)
-  aidsdeaths_noart <- array(sapply(dpsub("<AIDSDeathsNoART MV2>", c(4:20, 22:38), timedat.idx), as.numeric),
-                            lengths(dn5), dn5)
+  aidsdeaths_art5 <- array(sapply(dpsub("<AIDSDeathsART MV2>", c(4:20, 22:38), timedat.idx), as.numeric),
+                           lengths(dn5), dn5)
+  aidsdeaths_noart5 <- array(sapply(dpsub("<AIDSDeathsNoART MV2>", c(4:20, 22:38), timedat.idx), as.numeric),
+                             lengths(dn5), dn5)
 
 
   specres <- list("totpop.m" = totpop.m,
@@ -257,8 +257,8 @@ read_hivproj_output <- function(pjnz, single.age=TRUE){
                   "aidsdeaths.m" = aidsdeaths.m,
                   "aidsdeaths.f" = aidsdeaths.f,
                   aidsdeaths5 = aidsdeaths5,
-                  aidsdeaths_art = aidsdeaths_art,
-                  aidsdeaths_noart = aidsdeaths_noart)
+                  aidsdeaths_art5 = aidsdeaths_art5,
+                  aidsdeaths_noart5 = aidsdeaths_noart5)
 
   if(single.age){
 
@@ -305,12 +305,53 @@ read_hivproj_output <- function(pjnz, single.age=TRUE){
     else
       infections <- NA
 
+
+    ## Rob Glaubius, 25 October 2024
+    ## New tag <NonAIDSExcessDeathsSingleAge MV> stores estimated non-AIDS
+    ## excess deaths by:
+    ## * year (columns)
+    ## * sex (slowest-changing; both, male, female),
+    ## * single age (next slowest; all, 0…80+), and 
+    ## * ART status (fastest-moving; off/on).
+    ## 
+    ## Align output arrays similarly to other single-age outputs: age x sex x year
+    ## With 'male', and 'female', but not both sexes
+
+    ## 3 skips header; 82*2 skips 'both sex' block
+    noart_excess_mf_idx <- 3 + 82*2 + c(1:81, 83:163)*2
+    if (exists_dptag("<NonAIDSExcessDeathsSingleAge MV>")) {
+      nonaids_excess_deaths_noart <- sapply(dpsub("<NonAIDSExcessDeathsSingleAge MV>", noart_excess_mf_idx, timedat.idx), as.numeric)
+      nonaids_excess_deaths_art <- sapply(dpsub("<NonAIDSExcessDeathsSingleAge MV>", noart_excess_mf_idx + 1, timedat.idx), as.numeric)
+    } else {
+      nonaids_excess_deaths_noart <- rep(0.0, prod(lengths(dn1)))
+      nonaids_excess_deaths_art <- rep(0.0, prod(lengths(dn1)))
+    }
+
+
+    aidsdeath_row_idx <- 4 + 82 + c(0:80, 82+0:80)
+    if (exists_dptag("<AIDSDeathsARTSingleAge MV>")) {
+      aidsdeaths_art1 <- sapply(dpsub("<AIDSDeathsARTSingleAge MV>", aidsdeath_row_idx, timedat.idx), as.numeric)
+    } else {
+      aidsdeaths_art1 <- rep(0.0, prod(lengths(dn1)))
+    }
+
+    if (exists_dptag("<AIDSDeathsNoARTSingleAge MV>")) {
+      aidsdeaths_noart1 <- sapply(dpsub("<AIDSDeathsNoARTSingleAge MV>", aidsdeath_row_idx, timedat.idx), as.numeric)
+    } else {
+      aidsdeaths_noart1 <- rep(0.0, prod(lengths(dn1)))
+    }
+
     specres$totpop <- array(totpop, lengths(dn1), dn1)
     specres$hivpop <- array(hivpop, lengths(dn1), dn1)
     specres$artpop <- array(artpop, lengths(dn1), dn1)
     specres$natdeaths <- array(natdeaths, lengths(dn1), dn1)
     specres$hivdeaths <- array(hivdeaths, lengths(dn1), dn1)
     specres$infections <- array(infections, lengths(dn1), dn1)
+    specres$nonaids_excess_deaths_noart <- array(nonaids_excess_deaths_noart, lengths(dn1), dn1)
+    specres$nonaids_excess_deaths_art <- array(nonaids_excess_deaths_art, lengths(dn1), dn1)
+    specres$aidsdeaths_art1 <- array(aidsdeaths_art1, lengths(dn1), dn1)
+    specres$aidsdeaths_noart1 <- array(aidsdeaths_noart1, lengths(dn1), dn1)
+
   }
 
   specres$births <- stats::setNames(as.numeric(dpsub("<Births MV>", 2, timedat.idx)), proj.years)
@@ -563,6 +604,35 @@ read_hivproj_param <- function(pjnz, use_ep5=FALSE){
     artmx_timerr["ART1YR", ] <- val[2, ]
   }
 
+  ## Non-AIDS excess mortality by CD4
+  ## * Added in Spectrum 6.37 beta 17
+  ## * Initiated to default 0.0; will update witgh values from .DP if tag <AdultNonAIDSExcessMort MV> exists
+  ##
+  ## Formatting note from Rob Glaubius:
+  ## New tag <AdultNonAIDSExcessMort MV> stores the new rates.
+  ## These are organized into four rows for
+  ##   1) men off ART,
+  ##   2) men on ART,
+  ##   3) women off ART,
+  ##   4) women on ART.
+  ## 
+  ## Each row is laid out left-to-right in the same as our other adult HIV-related mortality rates:
+  ## * 15-24: CD4>500, 350-500, …, <50
+  ## * 25-34: CD4>500, 350-500, …, <50
+  ## * 35-44: CD4>500, 350-500, …, <50
+  ## * 45-54: CD4>500, 350-500, …, <50
+  ##
+
+  cd4_nonaids_excess_mort <- array(0.0, c(DS, 4, NG), dimnames(cd4_mort))
+  art_nonaids_excess_mort <- array(0.0, c(DS, 4, NG), dimnames(cd4_mort))
+
+  if(exists_dptag("<AdultNonAIDSExcessMort MV>")) {
+    cd4_nonaids_excess_mort[,,"Male"] <- array(as.numeric(dpsub("<AdultNonAIDSExcessMort MV>", 2, 4:31)), c(DS, 4))
+    art_nonaids_excess_mort[,,"Male"] <- array(as.numeric(dpsub("<AdultNonAIDSExcessMort MV>", 3, 4:31)), c(DS, 4))
+    cd4_nonaids_excess_mort[,,"Female"] <- array(as.numeric(dpsub("<AdultNonAIDSExcessMort MV>", 4, 4:31)), c(DS, 4))
+    art_nonaids_excess_mort[,,"Female"] <- array(as.numeric(dpsub("<AdultNonAIDSExcessMort MV>", 5, 4:31)), c(DS, 4))    
+  }
+
   ## program parameters
   if(dp.vers %in% c("<General 3>", "<General5>")){
     art15plus_numperc <- sapply(dp[adult.artnumperc.tidx+3:4, timedat.idx], as.numeric)
@@ -760,6 +830,8 @@ read_hivproj_param <- function(pjnz, use_ep5=FALSE){
                 "cd4_mort" = cd4_mort,
                 "art_mort" = art_mort,
                 "artmx_timerr" = artmx_timerr,
+                cd4_nonaids_excess_mort = cd4_nonaids_excess_mort,
+                art_nonaids_excess_mort = art_nonaids_excess_mort,
                 "art15plus_numperc" = art15plus_numperc,
                 "art15plus_num" = art15plus_num,
                 "adult_artadj_factor" = adult_artadj_factor,

diff --git a/R/spectrum.R b/R/spectrum.R
@@ -180,6 +180,10 @@ create_spectrum_fixpar <- function(projp, demp, hiv_steps_per_year = 10L, proj_s
   fp$art_mort <- projp$art_mort[c(1, 2, rep(3, hTS - 2)),,projp.h.ag,]
   fp$artmx_timerr <- projp$artmx_timerr[c(1, 2, rep(3, hTS - 2)), ]
 
+  fp$cd4_nonaids_excess_mort <- projp$cd4_nonaids_excess_mort[,projp.h.ag,]
+  fp$art_nonaids_excess_mort <- array(0.0, dim(fp$art_mort))
+  fp$art_nonaids_excess_mort[] <- rep(projp$art_nonaids_excess_mort[,projp.h.ag,], each = hTS)
+
   frr_agecat <- as.integer(rownames(projp$fert_rat))
   frr_agecat[frr_agecat == 18] <- 17
   fert_rat.h.ag <- findInterval(AGE_START + cumsum(h.ag.span[h.fert.idx]) - h.ag.span[h.fert.idx], frr_agecat)