2025 update -- non-AIDS excess mortality

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: first90
-Version: 1.6.11
+Version: 1.7.0
 Title: The first90 model
 Description: Implements the Shiny90 model for estimating progress towards the UNAIDS "first 90" target for HIV awareness of status in sub-Saharan Africa.
 Authors@R:

NEWS.md

@@ -1,3 +1,37 @@
+# first90 1.7.0
+
+* Update output plots and tables to 2024.
+
+* Extend theta0 initial values for additional year random-walk parameters.
+
+* 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.
+
+  `simmod()` 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.
+
+* Fix bug in R version for end-year net migration in first year of ART eligibility
+
+# first90 1.6.12
+
+* Update plots and tables to end in 2023.
+
 # first90 1.6.11
 
 * Fix bug in R version of simmod() where new infections previously tested negative were

R/create_fp.R

@@ -131,6 +131,10 @@ create_fp <- function(projp,
   fp$art_mort <- projp$art_mort[,,projp.h.ag,]
   fp$artmx_timerr <- projp$artmx_timerr
 
+  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)

R/create_hts_param.R

@@ -48,6 +48,8 @@ create_hts_param <- function(theta, fp) {
                      2022
                    } else if (length(theta) == 49) {
                      2023
+                   } else if (length(theta) == 51) {
+                     2024
                    } else {
                      stop("Unexpected length of parameter vector.")
                    }

R/eppasm.R

@@ -48,6 +48,17 @@ simmod <- 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))
+
+
   hivpopdeaths <- array(0, c(hDS, hAG, NG, PROJ_YEARS))
   artpopdeaths <- array(0, c(hTS, hDS, hAG, NG, PROJ_YEARS))
 
@@ -207,12 +218,21 @@ simmod <- 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 <- cd4_mort_ts * hivpop[,,,i]
+      hivdeaths.ts <- cd4_mort_ts * hivpop[,,,i]  # HIV mortality, untreated
+      grad <- grad - hivdeaths.ts              
       hivdeaths_hAG.ts <- colSums(hivdeaths.ts)
+      aidsdeaths_noart[,,,i] <- aidsdeaths_noart[,,,i] + DT * hivdeaths.ts
+
+      ## Non-AIDS excess mortality
+      nonaids_excess.ts <- fp$cd4_nonaids_excess_mort * hivpop[,,,i]
+      grad <- grad - nonaids_excess.ts
+      nonaids_excess_hAG.ts <- DT * colSums(nonaids_excess.ts)
+      excessnonaidsdeaths_noart[,,,i] <- excessnonaidsdeaths_noart[,,,i] + DT * nonaids_excess.ts
+
 
       ## ---- Distributing new infections in disease model ----
       if(fp$eppmod == "directinfections_hts") {
@@ -260,6 +280,9 @@ simmod <- function(fp, VERSION = "C") {
         prop_tn_hivp[!is.finite(prop_tn_hivp)] <- 0.0
         grad_tn[ , , hivp.idx] <- grad_tn[ , , hivp.idx] - hivdeaths_hAG.ts * prop_tn_hivp
 
+        ## Remove non-AIDS excess deaths among tested negative pop
+        grad_tn[ , , hivp.idx] <- grad_tn[ , , hivp.idx] - excess_nonaids_hAG.ts * prop_tn_hivp
+
         undiagnosed_i <- (hivpop[,,,i] - diagnpop[,,,i])
         prop_testneg <- testnegpop[ , , hivp.idx, i] / colSums(undiagnosed_i)
         prop_testneg[is.na(prop_testneg) | prop_testneg > 1 | prop_testneg < 0] <- 0
@@ -283,6 +306,7 @@ simmod <- function(fp, VERSION = "C") {
         grad_diagn[-1,,] <- grad_diagn[-1,,] + fp$cd4_prog * diagnpop[-hDS,,,i]      # add cd4 stage progression (untreated)
 
         grad_diagn <- grad_diagn - fp$cd4_mort * diagnpop[,,,i]                      # HIV mortality, untreated
+        grad_diagn <- grad_diagn - fp$cd4_nonaids_excess_mort * diagnpop[,,,i]       # Non-AIDS excess mortality
 
         diagnoses[,,,i] <- diagnoses[,,,i] + DT * (diagn_naive + diagn_testneg)
 
@@ -291,7 +315,7 @@ simmod <- function(fp, VERSION = "C") {
       }
 
       hivpop[,,,i] <- hivpop[,,,i] + DT*grad
-      hivpopdeaths[,,, i] <- hivpopdeaths[,,, i] + DT * hivdeaths.ts
+      hivpopdeaths[,,, i] <- hivpopdeaths[,,, i] + DT * (hivdeaths.ts + nonaids_excess.ts)
 
 
       ## ART population
@@ -308,8 +332,15 @@ simmod <- function(fp, VERSION = "C") {
 
         hivdeaths_hAG.ts <- hivdeaths_hAG.ts + colSums(artdeaths.ts,,2)
 
+        nonaids_excess_onart.ts <- fp$art_nonaids_excess_mort * artpop[,,,,i]
+        gradART <- gradART - nonaids_excess_onart.ts
+        nonaids_excess_hAG.ts <- nonaids_excess_hAG.ts +
+          DT * colSums(nonaids_excess_onart.ts,,2)
+        excessnonaidsdeaths_art[,,,,i] <- excessnonaidsdeaths_art[,,,,i] +
+          DT * nonaids_excess_onart.ts
+
         artpop[,,,, i] <- artpop[,,,, i] + DT * gradART
-        artpopdeaths[,,,, i] <- artpopdeaths[,,,, i] + DT * artdeaths.ts
+        artpopdeaths[,,,, i] <- artpopdeaths[,,,, i] + DT * (artdeaths.ts + nonaids_excess_onart.ts)
 
         ## ART dropout
         ## remove proportion from all adult ART groups back to untreated pop
@@ -496,6 +527,10 @@ simmod <- function(fp, VERSION = "C") {
       pop[,,hivp.idx,i] <- pop[,,hivp.idx,i] - hivdeaths_p.ts
       hivdeaths[,,i] <- hivdeaths[,,i] + hivdeaths_p.ts
 
+      nonaids_excess_p.ts <- apply(nonaids_excess_hAG.ts, 2, rep, h.ag.span) * apply(pop[,,hivp.idx,i], 2, calc.agdist)  # Non-AIDS excess deaths by single-year age
+      pop[,,hivp.idx,i] <- pop[,,hivp.idx,i] - nonaids_excess_p.ts
+      excessnonaidsdeaths[,,i] <- excessnonaidsdeaths[,,i] + nonaids_excess_p.ts      
+
     }
     # ---- End Disease Model ----
 
@@ -543,21 +578,21 @@ simmod <- function(fp, VERSION = "C") {
       hiv.mr.prob <- apply(mr.prob * pop[,,hivp.idx,i], 2, ctapply, ag.idx, sum) /  apply(pop[,,hivp.idx,i], 2, ctapply, ag.idx, sum)
       hiv.mr.prob[is.nan(hiv.mr.prob)] <- 0
 
-      if(i > fp$t_hts_start) {
+      if(i >= fp$t_hts_start) {
         hivn.mr.prob <- apply(mr.prob * pop[,,hivn.idx,i], 2, ctapply, ag.idx, sum) /  apply(pop[,,hivn.idx,i], 2, ctapply, ag.idx, sum)
         hivn.mr.prob[is.nan(hivn.mr.prob)] <- 0
       }
 
       pop[,,,i] <- sweep(pop[,,,i], 1:2, mr.prob, "*")
 
       hivpop[,,,i] <- sweep(hivpop[,,,i], 2:3, hiv.mr.prob, "*")
-      if(i > fp$t_hts_start) {
+      if(i >= fp$t_hts_start) {
         testnegpop[,, hivn.idx,i] <- testnegpop[,,hivn.idx,i] * hivn.mr.prob
         testnegpop[,, hivp.idx,i] <- testnegpop[,,hivp.idx,i] * hiv.mr.prob
 
         diagnpop[,,,i] <- sweep(diagnpop[,,,i], 2:3, hiv.mr.prob, "*")
       }
-      if(i > fp$tARTstart)
+      if(i >= fp$tARTstart)
         artpop[,,,,i] <- sweep(artpop[,,,,i], 3:4, hiv.mr.prob, "*")
     }
 
@@ -622,6 +657,17 @@ simmod <- function(fp, VERSION = "C") {
   attr(pop, "hivpopdeaths") <- hivpopdeaths
   attr(pop, "artpopdeaths") <- artpopdeaths
 
+  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, "hivtests") <- hivtests
   attr(pop, "diagnoses") <- diagnoses
   attr(pop, "late_diagnoses") <- late_diagnoses

R/extract-pjnz.R

@@ -101,6 +101,9 @@ extract_pjnz <- function(pjnz = NULL, dp_file= NULL, pjn_file = NULL) {
   v$cd4_mort <- get_dp_cd4_mort(dp)
   v$art_mort <- get_dp_art_mort(dp)
   v$artmx_timerr <- get_dp_artmx_timerr(dp, proj_years)
+
+  ## # Excess non-AIDS mortality
+  v <- c(v, get_dp_nonaids_excessmort(dp))
 
   ## # ART programme data
   v$art15plus_numperc <- array(as.numeric(unlist(dpsub(dp, "<HAARTBySexPerNum MV>", 4:5, col_idx))), lengths(dn), dn)
@@ -607,6 +610,45 @@ get_dp_art_mort <- function(dp) {
   art_mort
 }
 
+get_dp_nonaids_excessmort <- function(dp) {
+
+  ## 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
+  ##
+
+  dn <- list(cd4stage=1:DS,
+             agecat=c("15-24", "25-34", "35-44", "45+"),
+             sex=c("male", "female"))
+
+  val <- list()
+  val$cd4_nonaids_excess_mort <- array(0.0, c(DS, 4, NG), dn)
+  val$art_nonaids_excess_mort <- array(0.0, c(DS, 4, NG), dn)
+
+  if(exists_dptag(dp, "<AdultNonAIDSExcessMort MV>")) {
+    val$cd4_nonaids_excess_mort[,,"male"] <- array(as.numeric(dpsub(dp, "<AdultNonAIDSExcessMort MV>", 2, 4:31)), c(DS, 4))
+    val$art_nonaids_excess_mort[,,"male"] <- array(as.numeric(dpsub(dp, "<AdultNonAIDSExcessMort MV>", 3, 4:31)), c(DS, 4))
+    val$cd4_nonaids_excess_mort[,,"female"] <- array(as.numeric(dpsub(dp, "<AdultNonAIDSExcessMort MV>", 4, 4:31)), c(DS, 4))
+    val$art_nonaids_excess_mort[,,"female"] <- array(as.numeric(dpsub(dp, "<AdultNonAIDSExcessMort MV>", 5, 4:31)), c(DS, 4))    
+  }
+
+  val
+}
+
 get_dp_age14hivpop <- function(dp, proj_years) {
 
   PAED_DS <- 6 # number of paediatric stages of infection

R/inputs.R

@@ -37,7 +37,7 @@ select_prgmdata <- function(prgm_dat, cnt, age_group) {
     ## * year vector needs to be extended to output results to current year
 
     prg_dat <- data.frame(country = cnt,
-                          year = 2010:2022,
+                          year = 2010:2024,
                           agegr = '15-99', sex = 'both',
                           tot = NA, totpos = NA,
                           vct = NA, vctpos = NA, anc = NA, ancpos = NA)

R/likelihood.R

@@ -64,7 +64,7 @@ ll_prgdat <- function(mod, fp, dat) {
 ## -- UPDATE HERE --
 ## * max_year = <current_year> incremented each year
 
-art_constraint_penalty <- function(mod, fp, max_year = 2022) {
+art_constraint_penalty <- function(mod, fp, max_year = 2024) {
   ind_year <- c(2000:max_year) - fp$ss$proj_start + 1L
   tot_late <- apply(attr(mod, "late_diagnoses")[,,, ind_year], 4, sum)
   tot_untreated_pop <- apply(attr(mod, "hivpop")[,,, ind_year], 4, sum)
@@ -74,7 +74,7 @@ art_constraint_penalty <- function(mod, fp, max_year = 2022) {
   return(penalty)
 }
 # Include this in ll_hts if you want to incorporate the likelihood constraint on ART.
-  # val_art_penalty <- art_constraint_penalty(mod, fp, max_year = 2022)
+  # val_art_penalty <- art_constraint_penalty(mod, fp, max_year = 2024)
   # val <- val1 + val2 + val3 + val_prior + val_art_penalty
 
 # Function to prepare the data for input in the likelihood function.
@@ -142,7 +142,7 @@ lprior_hts <- function(theta, mod, fp) {
 
   ## -- UPDATE HERE --
   ## * Extend knots by 1 year to current year
-  knots <- 2000:2023 - fp$ss$proj_start + 1L
+  knots <- 2000:2024 - fp$ss$proj_start + 1L
   ## -- UPDATE ABOVE --
 
   n_k1 <- length(knots)

R/plot_functions.R

@@ -19,7 +19,7 @@ get_pjnz_summary_data <- function(fp) {
 #'
 ## -- UPDATE HERE --
 ## * update yr_pred to current year
-plot_pjnz_prv <- function(pjnz_summary, yr_pred = 2022) {
+plot_pjnz_prv <- function(pjnz_summary, yr_pred = 2024) {
   pjnz_summary <- stats::na.omit(data.frame(year = pjnz_summary[["year"]], prv = pjnz_summary[["prevalence"]]*100))
   pjnz_summary$year <- pjnz_summary$year + 0.5
   plot(pjnz_summary$prv ~ pjnz_summary$year,
@@ -34,7 +34,7 @@ plot_pjnz_prv <- function(pjnz_summary, yr_pred = 2022) {
 #' @export
 ## -- UPDATE HERE --
 ## * update yr_pred to current year
-plot_pjnz_inc <- function(pjnz_summary, yr_pred = 2022) {
+plot_pjnz_inc <- function(pjnz_summary, yr_pred = 2024) {
   pjnz_summary <- stats::na.omit(data.frame(year = pjnz_summary[["year"]], inc = pjnz_summary[["incidence"]]*1000))
   pjnz_summary$year <- pjnz_summary$year + 0.5
   pjnz_summary <- subset(pjnz_summary, year >= 2000)
@@ -50,7 +50,7 @@ plot_pjnz_inc <- function(pjnz_summary, yr_pred = 2022) {
 #' @export
 ## -- UPDATE HERE --
 ## * update yr_pred to current year
-plot_pjnz_pop <- function(pjnz_summary, yr_pred = 2022) {
+plot_pjnz_pop <- function(pjnz_summary, yr_pred = 2024) {
   pjnz_summary <- stats::na.omit(data.frame(year = pjnz_summary[["year"]], pop = pjnz_summary[["pop"]]/1000))
   pjnz_summary$year <- pjnz_summary$year + 0.5
   plot(pjnz_summary$pop ~ pjnz_summary$year,
@@ -65,7 +65,7 @@ plot_pjnz_pop <- function(pjnz_summary, yr_pred = 2022) {
 #' @export
 ## -- UPDATE HERE --
 ## * update yr_pred to current year
-plot_pjnz_plhiv <- function(pjnz_summary, yr_pred = 2022) {
+plot_pjnz_plhiv <- function(pjnz_summary, yr_pred = 2024) {
   pjnz_summary <- stats::na.omit(data.frame(year = pjnz_summary[["year"]], plhiv = pjnz_summary[["plhiv"]]/1000))
   pjnz_summary$year <- pjnz_summary$year + 0.5
   plot(pjnz_summary$plhiv ~ pjnz_summary$year,
@@ -83,7 +83,7 @@ plot_pjnz_plhiv <- function(pjnz_summary, yr_pred = 2022) {
 #' @export
 ## -- UPDATE HERE --
 ## * update yr_pred to current year
-plot_pjnz <- function(fp, yr_pred = 2022) {
+plot_pjnz <- function(fp, yr_pred = 2024) {
   summary <- get_pjnz_summary_data(fp)
   graphics::par(mfrow=c(2,2))
   plot_pjnz_prv(summary, yr_pred)
@@ -126,7 +126,7 @@ combine_rows <- function(prgdat) {
 #' @export
 ## -- UPDATE HERE --
 ## * update yr_pred to current year
-plot_input_tot <- function(prgdat, fp, yr_pred = 2022) {
+plot_input_tot <- function(prgdat, fp, yr_pred = 2024) {
     start <- fp$ss$proj_start
     mod <- simmod(fp)
     pop <- apply(mod[1:35,,,], 4, FUN=sum)
@@ -157,7 +157,7 @@ plot_input_tot <- function(prgdat, fp, yr_pred = 2022) {
 #' @export
 ## -- UPDATE HERE --
 ## * update yr_pred to current year
-plot_input_totpos <- function(prgdat, fp, yr_pred = 2022) {
+plot_input_totpos <- function(prgdat, fp, yr_pred = 2024) {
     start <- fp$ss$proj_start
     mod <- simmod(fp)
     plhiv <- apply(attr(mod, "hivpop")[,1:8,,], 4, FUN = sum) +
@@ -201,7 +201,7 @@ plot_input_totpos <- function(prgdat, fp, yr_pred = 2022) {
 #' @export
 ## -- UPDATE HERE --
 ## * update yr_pred to current year
-plot_input_anctot <- function(prgdat, fp, yr_pred = 2022) {
+plot_input_anctot <- function(prgdat, fp, yr_pred = 2024) {
 
   if (sum(prgdat$anc, na.rm = TRUE) > 0) {
     prgdat <- subset(prgdat, sex != 'male')
@@ -217,7 +217,7 @@ plot_input_anctot <- function(prgdat, fp, yr_pred = 2022) {
 #' @export
 ## -- UPDATE HERE --
 ## * update yr_pred to current year
-plot_input_ancpos <- function(prgdat, fp, yr_pred = 2022) {
+plot_input_ancpos <- function(prgdat, fp, yr_pred = 2024) {
 
   if (sum(prgdat$ancpos, na.rm = TRUE) > 0) {
     prgdat <- subset(prgdat, sex != 'male')
@@ -235,7 +235,7 @@ plot_input_ancpos <- function(prgdat, fp, yr_pred = 2022) {
 #' @export
 ## -- UPDATE HERE --
 ## * update yr_pred to current year
-plot_inputdata <- function(prgm_dat, fp, yr_pred = 2022) {
+plot_inputdata <- function(prgm_dat, fp, yr_pred = 2024) {
   graphics::par(mfrow = c(2,2))
   plot_input_tot(prgm_dat, fp, yr_pred)
   plot_input_totpos(prgm_dat, fp, yr_pred)
@@ -330,7 +330,7 @@ optimized_par <- function(opt, param = NULL) {
                    'RR re-testing: PLHIV aware (not ART) 2010',
 ## -- UPDATE HERE --
 ## * update label to current year
-                   'RR re-testing: PLHIV aware (not ART) 2022',
+                   'RR re-testing: PLHIV aware (not ART) 2024',
                    'RR re-testing: PLHIV on ART (*RR not ART)',
                    'RR among 25-34 men',
                    'RR among 35+ men',