Merge remote-tracking branch 'origin/master' into excess-nonaids-mx

mrc-ide · Jun 25, 2024 · fcb3d00 · fcb3d00
2 parents 7f34ae0 + 52688b2
commit fcb3d00
Show file tree

Hide file tree

Showing 34 changed files with 832 additions and 554 deletions.
diff --git a/.travis.yml b/.travis.yml
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -11,6 +11,8 @@ Imports:
     epp (>= 0.4.1),
     fastmatch (>= 1.1),
     mvtnorm,
+    plyr,
+    readxl,
     reshape2,
     vroom,
     xml2

diff --git a/NEWS.md b/NEWS.md
@@ -2,6 +2,18 @@
 
 * Add excess non-AIDS mortality among PLHIV. New model parameters `cd4_nonaids_excess_mort` and `art_nonaids_excess_mort`.
 
+## eppasm 0.7.6
+
+* Update internal data country ISO3 list to contain St. Kitts & Nevis and Dominica
+
+## eppasm 0.7.5
+
+* Qualify all package names and add all required packages into Imports section.
+
+## 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.

diff --git a/R/aggregate-fits.R b/R/aggregate-fits.R
@@ -22,7 +22,7 @@ create_aggr_input <- function(inputlist){
   ancrtcens <- do.call(rbind, lapply(eppdlist, "[[", "ancrtcens"))
   if(!is.null(ancrtcens) && nrow(ancrtcens)){
     ancrtcens$x <- ancrtcens$prev * ancrtcens$n
-    ancrtcens <- aggregate(cbind(x,n) ~ year, ancrtcens, sum)
+    ancrtcens <- stats::aggregate(cbind(x,n) ~ year, ancrtcens, sum)
     ancrtcens$prev <- ancrtcens$x / ancrtcens$n
     ancrtcens <- ancrtcens[c("year", "prev", "n")]
   }

diff --git a/R/cluster-functions.R b/R/cluster-functions.R
@@ -1,5 +1,5 @@
 
-pick_maxlpd <- function(set){set <- set[!sapply(set, is.null)]; set[[which.max(sapply(set, function(x) tail(x$stat,1)[1]))]]}
+pick_maxlpd <- function(set){set <- set[!sapply(set, is.null)]; set[[which.max(sapply(set, function(x) utils::tail(x$stat,1)[1]))]]}
 get_imisiter <- function(x) nrow(x$stat)
-get_logmargpost <- function(x) tail(x$stat, 1)[1]
-get_nunique <- function(x) tail(x$stat, 1)[2]
+get_logmargpost <- function(x) utils::tail(x$stat, 1)[1]
+get_nunique <- function(x) utils::tail(x$stat, 1)[2]
diff --git a/R/eppasm.R b/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
   }

diff --git a/R/fit-model.R b/R/fit-model.R
@@ -48,7 +48,7 @@ prepare_spec_fit <- function(pjnz, proj.end=2016.5, popadjust = NULL, popupdate=
 
 
   ## output
-  val <- setNames(vector("list", length(eppd)), names(eppd))
+  val <- stats::setNames(vector("list", length(eppd)), names(eppd))
 
   set.list.attr <- function(obj, attrib, value.lst)
     mapply(function(set, value){ attributes(set)[[attrib]] <- value; set}, obj, value.lst)
@@ -109,7 +109,7 @@ tidy_hhs_data <- function(eppd){
   hhs
 }
 
-create_subpop_specfp <- function(projp, demp, eppd, epp_t0=setNames(rep(1975, length(eppd)), names(eppd)), ..., popadjust=TRUE, popupdate=TRUE, perc_urban=NULL){
+create_subpop_specfp <- function(projp, demp, eppd, epp_t0=stats::setNames(rep(1975, length(eppd)), names(eppd)), ..., popadjust=TRUE, popupdate=TRUE, perc_urban=NULL){
 
   country <- attr(eppd, "country")
   country_code <- attr(eppd, "country_code")
@@ -195,7 +195,7 @@ create_subpop_specfp <- function(projp, demp, eppd, epp_t0=setNames(rep(1975, le
   subpop.dist <- prop.table(sapply(demp.subpop, get15to49pop, 2010))
 
   if(nrow(subset(eppd[[1]]$hhs, used)) != 0){ # HH survey data available
-    hhsprev.means <- sapply(lapply(eppd, function(dat) na.omit(dat$hhs$prev[dat$hhs$used])), mean)
+    hhsprev.means <- sapply(lapply(eppd, function(dat) stats::na.omit(dat$hhs$prev[dat$hhs$used])), mean)
     art.dist <- prop.table(subpop.dist * hhsprev.means)
   } else {  ## no HH survey data
     ## Apportion ART according to relative average ANC prevalence in each subpopulation
@@ -244,7 +244,7 @@ prepare_national_fit <- function(pjnz, upd.path=NULL, proj.end=2013.5, hiv_steps
   epp.input <- epp::read_epp_input(pjnz)
 
   ## output
-  val <- setNames(vector("list", length(eppd)), names(eppd))
+  val <- stats::setNames(vector("list", length(eppd)), names(eppd))
   val <- list()
 
   ## aggregate census data across regions
@@ -253,7 +253,7 @@ prepare_national_fit <- function(pjnz, upd.path=NULL, proj.end=2013.5, hiv_steps
     ancrtcens <- subset(ancrtcens, !is.na(prev) & !is.na(n))
     if(nrow(ancrtcens)){
       ancrtcens$x <- ancrtcens$prev * ancrtcens$n
-      ancrtcens <- aggregate(cbind(x,n) ~ year, ancrtcens, sum)
+      ancrtcens <- stats::aggregate(cbind(x,n) ~ year, ancrtcens, sum)
       ancrtcens$prev <- ancrtcens$x / ancrtcens$n
       ancrtcens <- ancrtcens[c("year", "prev", "n")]
     }
@@ -284,9 +284,9 @@ fitmod <- function(obj, ..., epp=FALSE, B0 = 1e5, B = 1e4, B.re = 3000, number_k
   ## ... : updates to fixed parameters (fp) object to specify fitting options
 
   if(epp)
-    fp <- update(attr(obj, 'eppfp'), ...)
+    fp <- stats::update(attr(obj, 'eppfp'), ...)
   else
-    fp <- update(attr(obj, 'specfp'), ...)
+    fp <- stats::update(attr(obj, 'specfp'), ...)
 
 
   ## Prepare likelihood data
@@ -347,13 +347,13 @@ fitmod <- function(obj, ..., epp=FALSE, B0 = 1e5, B = 1e4, B.re = 3000, number_k
       lpost0 <- likelihood(X0, fp, likdat, log=TRUE) + prior(X0, fp, log=TRUE)
       opt_init <- X0[which.max(lpost0)[1],]
     }
-    opt <- optim(opt_init, optfn, fp=fp, likdat=likdat, method=opt_method, control=list(fnscale=-1, trace=4, maxit=opt_maxit, ndeps=rep(opt_diffstep, length(opt_init))))
+    opt <- stats::optim(opt_init, optfn, fp=fp, likdat=likdat, method=opt_method, control=list(fnscale=-1, trace=4, maxit=opt_maxit, ndeps=rep(opt_diffstep, length(opt_init))))
     opt$fp <- fp
     opt$likdat <- likdat
     opt$param <- fnCreateParam(opt$par, fp)
-    opt$mod <- simmod(update(fp, list=opt$param))
+    opt$mod <- simmod(stats::update(fp, list=opt$param))
     if(opthess){
-      opt$hessian <- optimHess(opt_init, optfn, fp=fp, likdat=likdat,
+      opt$hessian <- stats::optimHess(opt_init, optfn, fp=fp, likdat=likdat,
                                control=list(fnscale=-1,
                                             trace=4,
                                             maxit=1e3,
@@ -450,7 +450,7 @@ simfit.specfit <- function(fit,
     if(rwproj)
       fit <- rw_projection(fit)
 
-    fp_list <- lapply(fit$param, function(par) update(fit$fp, list=par))
+    fp_list <- lapply(fit$param, function(par) stats::update(fit$fp, list=par))
     mod.list <- lapply(fp_list, simmod)
 
   } else {
@@ -598,7 +598,7 @@ simfit.eppfit <- function(fit, rwproj=fit$fp$eppmod == "rspline", pregprev=TRUE)
     fit$param <- lapply(fit$param, function(par){par$rvec <- sim_rvec_rwproj(par$rvec, firstidx, lastidx, fit$fp$dt); par})
   }
 
-  fp.list <- lapply(fit$param, function(par) update(fit$fp, list=par))
+  fp.list <- lapply(fit$param, function(par) stats::update(fit$fp, list=par))
   mod.list <- lapply(fp.list, simmod)
 
   fit$rvec <- sapply(mod.list, attr, "rvec")
@@ -636,7 +636,7 @@ sim_mod_list <- function(fit, rwproj=fit$fp$eppmod == "rspline"){
     fit$param <- lapply(fit$param, function(par){par$rvec <- epp:::sim_rvec_rwproj(par$rvec, firstidx, lastidx, dt); par})
   }
 
-  fp.list <- lapply(fit$param, function(par) update(fit$fp, list=par))
+  fp.list <- lapply(fit$param, function(par) stats::update(fit$fp, list=par))
   mod.list <- lapply(fp.list, simmod)
 
   ## strip unneeded attributes to preserve memory