diff --git a/DESCRIPTION b/DESCRIPTION index fe30c69..0a9e674 100644 --- 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 = "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 b936e37..9f292ab 100644 --- 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 `""`. + +* `read_hivproj_output()` updated to read additiounal deaths outputs: + - Non-AIDS excess deaths among PLHIV from the tag `""`, + 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 index 13f3aae..8f941df 100644 --- 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 index a90f3da..6662dc4 100644 --- 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("", c(4:20, 22:38), timedat.idx), as.numeric), - lengths(dn5), dn5) - aidsdeaths_noart <- array(sapply(dpsub("", c(4:20, 22:38), timedat.idx), as.numeric), - lengths(dn5), dn5) + aidsdeaths_art5 <- array(sapply(dpsub("", c(4:20, 22:38), timedat.idx), as.numeric), + lengths(dn5), dn5) + aidsdeaths_noart5 <- array(sapply(dpsub("", 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 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("")) { + nonaids_excess_deaths_noart <- sapply(dpsub("", noart_excess_mf_idx, timedat.idx), as.numeric) + nonaids_excess_deaths_art <- sapply(dpsub("", 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("")) { + aidsdeaths_art1 <- sapply(dpsub("", aidsdeath_row_idx, timedat.idx), as.numeric) + } else { + aidsdeaths_art1 <- rep(0.0, prod(lengths(dn1))) + } + + if (exists_dptag("")) { + aidsdeaths_noart1 <- sapply(dpsub("", 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("", 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 exists + ## + ## Formatting note from Rob Glaubius: + ## New tag 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("")) { + cd4_nonaids_excess_mort[,,"Male"] <- array(as.numeric(dpsub("", 2, 4:31)), c(DS, 4)) + art_nonaids_excess_mort[,,"Male"] <- array(as.numeric(dpsub("", 3, 4:31)), c(DS, 4)) + cd4_nonaids_excess_mort[,,"Female"] <- array(as.numeric(dpsub("", 4, 4:31)), c(DS, 4)) + art_nonaids_excess_mort[,,"Female"] <- array(as.numeric(dpsub("", 5, 4:31)), c(DS, 4)) + } + ## program parameters if(dp.vers %in% c("", "")){ 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 index 9d1a0b7..ab7a7e6 100644 --- 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) diff --git a/src/eppasm.cpp b/src/eppasm.cpp index 61b1283..e77184e 100644 --- a/src/eppasm.cpp +++ b/src/eppasm.cpp @@ -130,7 +130,9 @@ extern "C" { multi_array_ref cd4_initdist(REAL(getListElement(s_fp, "cd4_initdist")), extents[NG][hAG][hDS]); multi_array_ref cd4_prog(REAL(getListElement(s_fp, "cd4_prog")), extents[NG][hAG][hDS-1]); multi_array_ref cd4_mort(REAL(getListElement(s_fp, "cd4_mort")), extents[NG][hAG][hDS]); + multi_array_ref cd4_nonaids_excess_mort(REAL(getListElement(s_fp, "cd4_nonaids_excess_mort")), extents[NG][hAG][hDS]); multi_array_ref art_mort(REAL(getListElement(s_fp, "art_mort")), extents[NG][hAG][hDS][hTS]); + multi_array_ref art_nonaids_excess_mort(REAL(getListElement(s_fp, "art_nonaids_excess_mort")), extents[NG][hAG][hDS][hTS]); multi_array_ref artmx_timerr(REAL(getListElement(s_fp, "artmx_timerr")), extents[PROJ_YEARS][hTS]); // sub-fertility @@ -291,6 +293,16 @@ extern "C" { Rf_setAttrib(s_pop, Rf_install("natdeaths"), s_natdeaths); multi_array_ref natdeaths(REAL(s_natdeaths), extents[PROJ_YEARS][NG][pAG]); memset(REAL(s_natdeaths), 0, Rf_length(s_natdeaths)*sizeof(double)); + + SEXP s_excessnonaidsdeaths = PROTECT(Rf_allocVector(REALSXP, pAG * NG * PROJ_YEARS)); + SEXP s_excessnonaidsdeaths_dim = PROTECT(Rf_allocVector(INTSXP, 3)); + INTEGER(s_excessnonaidsdeaths_dim)[0] = pAG; + INTEGER(s_excessnonaidsdeaths_dim)[1] = NG; + INTEGER(s_excessnonaidsdeaths_dim)[2] = PROJ_YEARS; + Rf_setAttrib(s_excessnonaidsdeaths, R_DimSymbol, s_excessnonaidsdeaths_dim); + Rf_setAttrib(s_pop, Rf_install("excessnonaidsdeaths"), s_excessnonaidsdeaths); + multi_array_ref excessnonaidsdeaths(REAL(s_excessnonaidsdeaths), extents[PROJ_YEARS][NG][pAG]); + memset(REAL(s_excessnonaidsdeaths), 0, Rf_length(s_excessnonaidsdeaths)*sizeof(double)); SEXP s_aidsdeaths_noart = PROTECT(Rf_allocVector(REALSXP, hDS * hAG * NG * PROJ_YEARS)); SEXP s_aidsdeaths_noart_dim = PROTECT(Rf_allocVector(INTSXP, 4)); @@ -303,6 +315,28 @@ extern "C" { multi_array_ref aidsdeaths_noart(REAL(s_aidsdeaths_noart), extents[PROJ_YEARS][NG][hAG][hDS]); memset(REAL(s_aidsdeaths_noart), 0, Rf_length(s_aidsdeaths_noart)*sizeof(double)); + SEXP s_natdeaths_noart = PROTECT(Rf_allocVector(REALSXP, hDS * hAG * NG * PROJ_YEARS)); + SEXP s_natdeaths_noart_dim = PROTECT(Rf_allocVector(INTSXP, 4)); + INTEGER(s_natdeaths_noart_dim)[0] = hDS; + INTEGER(s_natdeaths_noart_dim)[1] = hAG; + INTEGER(s_natdeaths_noart_dim)[2] = NG; + INTEGER(s_natdeaths_noart_dim)[3] = PROJ_YEARS; + Rf_setAttrib(s_natdeaths_noart, R_DimSymbol, s_natdeaths_noart_dim); + Rf_setAttrib(s_pop, Rf_install("natdeaths_noart"), s_natdeaths_noart); + multi_array_ref natdeaths_noart(REAL(s_natdeaths_noart), extents[PROJ_YEARS][NG][hAG][hDS]); + memset(REAL(s_natdeaths_noart), 0, Rf_length(s_natdeaths_noart)*sizeof(double)); + + SEXP s_excessnonaidsdeaths_noart = PROTECT(Rf_allocVector(REALSXP, hDS * hAG * NG * PROJ_YEARS)); + SEXP s_excessnonaidsdeaths_noart_dim = PROTECT(Rf_allocVector(INTSXP, 4)); + INTEGER(s_excessnonaidsdeaths_noart_dim)[0] = hDS; + INTEGER(s_excessnonaidsdeaths_noart_dim)[1] = hAG; + INTEGER(s_excessnonaidsdeaths_noart_dim)[2] = NG; + INTEGER(s_excessnonaidsdeaths_noart_dim)[3] = PROJ_YEARS; + Rf_setAttrib(s_excessnonaidsdeaths_noart, R_DimSymbol, s_excessnonaidsdeaths_noart_dim); + Rf_setAttrib(s_pop, Rf_install("excessnonaidsdeaths_noart"), s_excessnonaidsdeaths_noart); + multi_array_ref excessnonaidsdeaths_noart(REAL(s_excessnonaidsdeaths_noart), extents[PROJ_YEARS][NG][hAG][hDS]); + memset(REAL(s_excessnonaidsdeaths_noart), 0, Rf_length(s_excessnonaidsdeaths_noart)*sizeof(double)); + SEXP s_aidsdeaths_art = PROTECT(Rf_allocVector(REALSXP, hTS * hDS * hAG * NG * PROJ_YEARS)); SEXP s_aidsdeaths_art_dim = PROTECT(Rf_allocVector(INTSXP, 5)); INTEGER(s_aidsdeaths_art_dim)[0] = hTS; @@ -315,6 +349,29 @@ extern "C" { multi_array_ref aidsdeaths_art(REAL(s_aidsdeaths_art), extents[PROJ_YEARS][NG][hAG][hDS][hTS]); memset(REAL(s_aidsdeaths_art), 0, Rf_length(s_aidsdeaths_art)*sizeof(double)); + SEXP s_natdeaths_art = PROTECT(allocVector(REALSXP, hTS * hDS * hAG * NG * PROJ_YEARS)); + SEXP s_natdeaths_art_dim = PROTECT(allocVector(INTSXP, 5)); + INTEGER(s_natdeaths_art_dim)[0] = hTS; + INTEGER(s_natdeaths_art_dim)[1] = hDS; + INTEGER(s_natdeaths_art_dim)[2] = hAG; + INTEGER(s_natdeaths_art_dim)[3] = NG; + INTEGER(s_natdeaths_art_dim)[4] = PROJ_YEARS; + setAttrib(s_natdeaths_art, R_DimSymbol, s_natdeaths_art_dim); + setAttrib(s_pop, install("natdeaths_art"), s_natdeaths_art); + multi_array_ref natdeaths_art(REAL(s_natdeaths_art), extents[PROJ_YEARS][NG][hAG][hDS][hTS]); + memset(REAL(s_natdeaths_art), 0, length(s_natdeaths_art)*sizeof(double)); + + SEXP s_excessnonaidsdeaths_art = PROTECT(allocVector(REALSXP, hTS * hDS * hAG * NG * PROJ_YEARS)); + SEXP s_excessnonaidsdeaths_art_dim = PROTECT(allocVector(INTSXP, 5)); + INTEGER(s_excessnonaidsdeaths_art_dim)[0] = hTS; + INTEGER(s_excessnonaidsdeaths_art_dim)[1] = hDS; + INTEGER(s_excessnonaidsdeaths_art_dim)[2] = hAG; + INTEGER(s_excessnonaidsdeaths_art_dim)[3] = NG; + INTEGER(s_excessnonaidsdeaths_art_dim)[4] = PROJ_YEARS; + setAttrib(s_excessnonaidsdeaths_art, R_DimSymbol, s_excessnonaidsdeaths_art_dim); + setAttrib(s_pop, install("excessnonaidsdeaths_art"), s_excessnonaidsdeaths_art); + multi_array_ref excessnonaidsdeaths_art(REAL(s_excessnonaidsdeaths_art), extents[PROJ_YEARS][NG][hAG][hDS][hTS]); + memset(REAL(s_excessnonaidsdeaths_art), 0, length(s_excessnonaidsdeaths_art)*sizeof(double)); SEXP s_popadjust = PROTECT(Rf_allocVector(REALSXP, pAG * NG * PROJ_YEARS)); SEXP s_popadjust_dim = PROTECT(Rf_allocVector(INTSXP, 3)); @@ -499,6 +556,7 @@ extern "C" { int a = 0; for(int ha = 0; ha < hAG; ha++){ double deathsmig_ha = 0, hivpop_ha = 0; + double deaths_ha = 0.0; for(int i = 0; i < hAG_SPAN[ha]; i++){ hivpop_ha += pop[t][HIVP][g][a]; @@ -508,6 +566,7 @@ extern "C" { double ndeaths_a = pop[t][HIVN][g][a] * qx; pop[t][HIVN][g][a] -= ndeaths_a; // survival HIV- population double hdeaths_a = pop[t][HIVP][g][a] * qx; + deaths_ha += hdeaths_a; deathsmig_ha -= hdeaths_a; pop[t][HIVP][g][a] -= hdeaths_a; // survival HIV+ population natdeaths[t][g][a] = ndeaths_a + hdeaths_a; @@ -524,12 +583,16 @@ extern "C" { } // migration and deaths for hivpop + double deathrate_ha = hivpop_ha > 0 ? deaths_ha / hivpop_ha : 0.0; double deathmigrate_ha = hivpop_ha > 0 ? deathsmig_ha / hivpop_ha : 0.0; for(int hm = 0; hm < hDS; hm++){ + natdeaths_noart[t][g][ha][hm] += hivpop[t][g][ha][hm] * deathrate_ha; hivpop[t][g][ha][hm] *= 1+deathmigrate_ha; if(t > t_ART_start) - for(int hu = 0; hu < hTS; hu++) + for(int hu = 0; hu < hTS; hu++) { + natdeaths_art[t][g][ha][hm][hu] += artpop[t][g][ha][hm][hu] * deathrate_ha; artpop[t][g][ha][hm][hu] *= 1+deathmigrate_ha; + } } // loop over hm } // loop over ha } // loop over g @@ -593,6 +656,9 @@ extern "C" { double hivdeaths_ha[NG][hAG]; memset(hivdeaths_ha, 0, sizeof(double)*NG*hAG); + double nonaids_excess_ha[NG][hAG]; + memset(nonaids_excess_ha, 0, sizeof(double)*NG*hAG); + // untreated population // disease progression and mortality @@ -610,10 +676,16 @@ extern "C" { 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]; - hivdeaths_ha[g][ha] += DT*deaths; - aidsdeaths_noart[t][g][ha][hm] += DT*deaths; - grad[g][ha][hm] = -deaths; + double aids_deaths = cd4mx_scale * cd4_mort[g][ha][hm] * hivpop[t][g][ha][hm]; + hivdeaths_ha[g][ha] += DT * aids_deaths; + aidsdeaths_noart[t][g][ha][hm] += DT * aids_deaths; + + double excess_nonaids_deaths = cd4_nonaids_excess_mort[g][ha][hm] * hivpop[t][g][ha][hm]; + nonaids_excess_ha[g][ha] += DT * excess_nonaids_deaths; + excessnonaidsdeaths_noart[t][g][ha][hm] += DT * excess_nonaids_deaths; + + + grad[g][ha][hm] = -(aids_deaths + excess_nonaids_deaths); } for(int hm = 1; hm < hDS; hm++){ grad[g][ha][hm-1] -= cd4_prog[g][ha][hm-1] * hivpop[t][g][ha][hm-1]; @@ -696,10 +768,15 @@ extern "C" { for(int hm = everARTelig_idx; hm < hDS; hm++){ for(int hu = 0; hu < hTS; hu++){ - double deaths = art_mort[g][ha][hm][hu] * artmx_timerr[t][hu] * artpop[t][g][ha][hm][hu]; - hivdeaths_ha[g][ha] += DT*deaths; - aidsdeaths_art[t][g][ha][hm][hu] += DT*deaths; - gradART[g][ha][hm][hu] = -deaths; + double aids_deaths = art_mort[g][ha][hm][hu] * artmx_timerr[t][hu] * artpop[t][g][ha][hm][hu]; + double nonaids_deaths = art_nonaids_excess_mort[g][ha][hm][hu] * artpop[t][g][ha][hm][hu]; + hivdeaths_ha[g][ha] += DT * aids_deaths; + aidsdeaths_art[t][g][ha][hm][hu] += DT * aids_deaths; + + nonaids_excess_ha[g][ha] += DT * nonaids_deaths; + excessnonaidsdeaths_art[t][g][ha][hm][hu] += DT * nonaids_deaths; + + gradART[g][ha][hm][hu] = -(aids_deaths + nonaids_deaths); } for(int hu = 0; hu < (hTS - 1); hu++) { @@ -969,9 +1046,11 @@ extern "C" { for(int ha = 0; ha < hAG; ha++){ if(hivpop_ha[ha] > 0){ double hivqx_ha = hivdeaths_ha[g][ha] / hivpop_ha[ha]; + double nonaids_excess_qx_ha = nonaids_excess_ha[g][ha] / hivpop_ha[ha]; for(int i = 0; i < hAG_SPAN[ha]; i++){ hivdeaths[t][g][a] += pop[t][HIVP][g][a] * hivqx_ha; - pop[t][HIVP][g][a] *= (1.0-hivqx_ha); + excessnonaidsdeaths[t][g][a] += pop[t][HIVP][g][a] * nonaids_excess_qx_ha; + pop[t][HIVP][g][a] *= (1.0 - hivqx_ha - nonaids_excess_qx_ha); a++; } } else { @@ -1127,7 +1206,7 @@ extern "C" { incid15to49[t] /= incid15to49_denom; } - UNPROTECT(28); + UNPROTECT(38); return s_pop; } }