Merge pull request #41 from mrc-ide/art-dropout-cd4-recovery

ART dropout cd4 recovery

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: eppasm
 Title: Age-structured EPP Model for HIV Epidemic Estimates
-Version: 0.7.3
+Version: 0.7.4
 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),

NEWS.md

@@ -1,3 +1,7 @@
+## eppasm 0.7.4
+
+* Implement recovery to next higher CD4 category following ART interruption for those on ART greater than one year.
+
 ## eppasm 0.7.3
 
 * Bug fix: account for end-year net migration in the ART population in the first year of ART start.

R/eppasm.R

@@ -1,15 +1,15 @@
 
 #' @useDynLib eppasm eppasmC
 #' @export
-simmod.specfp <- function(fp, VERSION="C", ...){
+simmod.specfp <- function(fp, VERSION="C", ...) {
 
   if(!exists("popadjust", where=fp))
     fp$popadjust <- FALSE
 
   if(!exists("incidmod", where=fp))
     fp$incidmod <- "eppspectrum"
 
-  if(VERSION != "R"){
+  if(VERSION != "R") {
 
     ## eppmod codes:
     ## 0: r-spline
@@ -29,7 +29,7 @@ simmod.specfp <- function(fp, VERSION="C", ...){
     return(mod)
   }
 
-##################################################################################
+  ##################################################################################
 
   if(requireNamespace("fastmatch", quietly = TRUE))
     ctapply <- fastmatch::ctapply
@@ -62,7 +62,7 @@ simmod.specfp <- function(fp, VERSION="C", ...){
 
   popadj.prob <- array(0, c(pAG, NG, PROJ_YEARS))
 
-  if(fp$eppmod != "directincid_ann"){
+  if(fp$eppmod != "directincid_ann") {
     ## outputs by timestep
     incrate15to49.ts.out <- rep(NA, length(fp$rvec))
     rvec <- if(fp$eppmod == "rtrend") rep(NA, length(fp$proj.steps)) else fp$rvec
@@ -90,7 +90,7 @@ simmod.specfp <- function(fp, VERSION="C", ...){
     ## Add lagged births into youngest age group
     entrant_prev <- fp$entrantprev[,i]
 
-    if(exists("popadjust", where=fp) & fp$popadjust){
+    if(exists("popadjust", where=fp) & fp$popadjust) {
       hivn_entrants <- fp$entrantpop[,i-1]*(1-entrant_prev)
       hivp_entrants <- fp$entrantpop[,i-1]*entrant_prev
     } else {
@@ -112,7 +112,7 @@ simmod.specfp <- function(fp, VERSION="C", ...){
     hivpop[,-1,,i] <- hivpop[,-1,,i] + sweep(hivpop[,-hAG,,i-1], 2:3, hiv.ag.prob[-hAG,], "*")
     hivpop[,1,,i] <- hivpop[,1,,i] + sweep(fp$paedsurv_cd4dist[,,i], 2, hivp_entrants * (1-fp$entrantartcov[,i]), "*")
 
-    if(i > fp$tARTstart){
+    if(i > fp$tARTstart) {
       artpop[,,,,i] <- artpop[,,,,i-1]
       artpop[,,-hAG,,i] <- artpop[,,-hAG,,i] - sweep(artpop[,,-hAG,,i-1], 3:4, hiv.ag.prob[-hAG,], "*")
       artpop[,,-1,,i] <- artpop[,,-1,,i] + sweep(artpop[,,-hAG,,i-1], 3:4, hiv.ag.prob[-hAG,], "*")
@@ -213,7 +213,7 @@ simmod.specfp <- function(fp, VERSION="C", ...){
       grad[-hDS,,] <- grad[-hDS,,] - fp$cd4_prog * hivpop[-hDS,,,i]  # remove cd4 stage progression (untreated)
       grad[-1,,] <- grad[-1,,] + fp$cd4_prog * hivpop[-hDS,,,i]      # add cd4 stage progression (untreated)
 
-      if(fp$scale_cd4_mort == 1){
+      if(fp$scale_cd4_mort == 1) {
         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
@@ -244,7 +244,18 @@ simmod.specfp <- function(fp, VERSION="C", ...){
 
         ## ART dropout
         ## remove proportion from all adult ART groups back to untreated pop
-        grad <- grad + fp$art_dropout[i]*colSums(artpop[,,,,i])
+        art_dropout_ii <- fp$art_dropout[i]*colSums(artpop[1:2,,,,i])
+        if (fp$art_dropout_recover_cd4) {
+          art_dropout_ii[1,,] <- art_dropout_ii[1,,] +
+            fp$art_dropout[i] * artpop[3:fp$ss$hTS,1,,,i]
+          art_dropout_ii[-fp$ss$hDS,,] <- art_dropout_ii[-fp$ss$hDS,,] +
+            fp$art_dropout[i] * artpop[3:fp$ss$hTS,-1,,,i]
+        } else {
+          art_dropout_ii <- art_dropout_ii +
+            fp$art_dropout[i] * artpop[3:fp$ss$hTS,,,,i]
+        }
+
+        grad <- grad + art_dropout_ii
         gradART <- gradART - fp$art_dropout[i]*artpop[,,,,i]
 
         ## calculate number eligible for ART
@@ -255,7 +266,7 @@ simmod.specfp <- function(fp, VERSION="C", ...){
         art15plus.elig <- sweep(hivpop[,h.age15plus.idx,,i], 1, artcd4_percelig, "*")
 
         ## calculate pregnant women
-        if(fp$pw_artelig[i]){
+        if(fp$pw_artelig[i]) {
           births.dist <- sweep(fp$frr_cd4[,,i] * hivpop[,h.fert.idx,f.idx,i], 2,
                                births.by.h.age / (ctapply(pop[p.fert.idx, f.idx, hivn.idx, i], ag.idx[p.fert.idx], sum) + colSums(fp$frr_cd4[,,i] * hivpop[,h.fert.idx,f.idx,i]) + colSums(fp$frr_art[,,,i] * artpop[ ,,h.fert.idx,f.idx,i],,2)), "*")
           if(fp$artcd4elig_idx[i] > 1)
@@ -266,14 +277,14 @@ simmod.specfp <- function(fp, VERSION="C", ...){
 
         artpop_curr_g <- colSums(artpop[,,h.age15plus.idx,,i],,3) + DT*colSums(gradART[,,h.age15plus.idx,],,3)
         artnum.ii <- c(0,0) # number on ART this ts
-        if (fp$projection_period == "midyear" && DT*ii < 0.5){
+        if (fp$projection_period == "midyear" && DT*ii < 0.5) {
           for(g in 1:2){
-            if(!any(fp$art15plus_isperc[g,i-2:1])){  # both number
+            if(!any(fp$art15plus_isperc[g,i-2:1])) {  # both number
               artnum.ii[g] <- c(fp$art15plus_num[g,i-2:1] %*% c(1-(DT*ii+0.5), DT*ii+0.5))
-            } else if(all(fp$art15plus_isperc[g,i-2:1])){  # both percentage
+            } else if(all(fp$art15plus_isperc[g,i-2:1])) {  # both percentage
               artcov.ii <- c(fp$art15plus_num[g,i-2:1] %*% c(1-(DT*ii+0.5), DT*ii+0.5))
               artnum.ii[g] <- artcov.ii * (sum(art15plus.elig[,,g]) + artpop_curr_g[g])
-            } else if(!fp$art15plus_isperc[g,i-2] & fp$art15plus_isperc[g,i-1]){ # transition number to percentage
+            } else if(!fp$art15plus_isperc[g,i-2] & fp$art15plus_isperc[g,i-1]) { # transition number to percentage
               curr_coverage <- artpop_curr_g[g] / (sum(art15plus.elig[,,g]) + artpop_curr_g[g])
               artcov.ii <- curr_coverage + (fp$art15plus_num[g,i-1] - curr_coverage) * DT/(0.5-DT*(ii-1))
               artnum.ii[g] <- artcov.ii * (sum(art15plus.elig[,,g]) + artpop_curr_g[g])
@@ -286,12 +297,12 @@ simmod.specfp <- function(fp, VERSION="C", ...){
               art_interp_w <- art_interp_w - 0.5
             }
 
-            if(!any(fp$art15plus_isperc[g,i-1:0])){  # both number
+            if(!any(fp$art15plus_isperc[g,i-1:0])) {  # both number
               artnum.ii[g] <- c(fp$art15plus_num[g,i-1:0] %*% c(1-art_interp_w, art_interp_w))
             } else if(all(fp$art15plus_isperc[g,i-1:0])) {  # both percentage
               artcov.ii <- c(fp$art15plus_num[g,i-1:0] %*% c(1-art_interp_w, art_interp_w))
               artnum.ii[g] <- artcov.ii * (sum(art15plus.elig[,,g]) + artpop_curr_g[g])
-            } else if(!fp$art15plus_isperc[g,i-1] & fp$art15plus_isperc[g,i]){  # transition number to percentage
+            } else if(!fp$art15plus_isperc[g,i-1] & fp$art15plus_isperc[g,i]) {  # transition number to percentage
               curr_coverage <- artpop_curr_g[g] / (sum(art15plus.elig[,,g]) + artpop_curr_g[g])
               artcov.ii <- curr_coverage + (fp$art15plus_num[g,i] - curr_coverage) * DT/(1.0 - art_interp_w + DT)              
               artnum.ii[g] <- artcov.ii * (sum(art15plus.elig[,,g]) + artpop_curr_g[g])
@@ -303,9 +314,9 @@ simmod.specfp <- function(fp, VERSION="C", ...){
         art15plus.inits <- pmax(artnum.ii - artpop_curr_g, 0)
 
         ## calculate ART initiation distribution
-        if(!fp$med_cd4init_input[i]){
+        if(!fp$med_cd4init_input[i]) {
 
-          if(fp$art_alloc_method == 4L){ ## by lowest CD4
+          if(fp$art_alloc_method == 4L) { ## by lowest CD4
 
             ## Calculate proportion to be initiated in each CD4 category
             artinit <- array(0, dim(art15plus.elig))
@@ -383,7 +394,7 @@ simmod.specfp <- function(fp, VERSION="C", ...){
     ## incrate15to49.i <- (fp$prev15to49[i] - prev.i)/(1-prev.i)
 
     ## Direct incidence input
-    if(fp$eppmod == "directincid_ann"){
+    if(fp$eppmod == "directincid_ann") {
       agesex.inc <- sweep(fp$incrr_age[,,i], 2, sexinc/(colSums(pop[p.incidpop.idx,,hivn.idx,i] * fp$incrr_age[p.incidpop.idx,,i])/colSums(pop[p.incidpop.idx,,hivn.idx,i-1])), "*")
       infections[,,i] <- agesex.inc * pop[,,hivn.idx,i]
       pop[,,hivn.idx,i] <- pop[,,hivn.idx,i] - infections[,,i]
@@ -409,7 +420,7 @@ simmod.specfp <- function(fp, VERSION="C", ...){
 
 
     ## adjust population to match target population size
-    if(exists("popadjust", where=fp) & fp$popadjust){
+    if(exists("popadjust", where=fp) & fp$popadjust) {
       popadj.prob[,,i] <- fp$targetpop[,,i] / rowSums(pop[,,,i],,2)
       hiv.popadj.prob <- apply(popadj.prob[,,i] * pop[,,2,i], 2, ctapply, ag.idx, sum) /  apply(pop[,,2,i], 2, ctapply, ag.idx, sum)
       hiv.popadj.prob[is.nan(hiv.popadj.prob)] <- 0
@@ -456,7 +467,7 @@ simmod.specfp <- function(fp, VERSION="C", ...){
 
   attr(pop, "pregprevlag") <- pregprevlag
 
-  if(fp$eppmod != "directincid_ann"){
+  if(fp$eppmod != "directincid_ann") {
     attr(pop, "incrate15to49_ts") <- incrate15to49.ts.out
     attr(pop, "prev15to49_ts") <- prev15to49.ts.out
   }