Merge pull request #31 from statnet/prep-sti

Merging current CID model into master

DESCRIPTION

@@ -1,11 +1,12 @@
 Package: EpiModelHIV
-Version: 1.0.0
-Date: 2016-06-25
+Version: 1.5.0
+Date: 2017-05-04
 Type: Package
-Title: Epidemic Modeling of HIV Transmission among MSM and Heterosexual Populations
+Title: Network-Based Epidemic Modeling of HIV Transmission among MSM and Heterosexual Populations
 Authors@R: c(person("Samuel M.", "Jenness", role = c("cre", "aut"), email = "samuel.m.jenness@emory.edu"),
              person("Steven M.", "Goodreau", role="aut", email="goodeau@uw.edu"),
-             person("Emily", "Beylerian", role = "ctb", email = "ebey@uw.edu"))
+             person("Emily", "Beylerian", role = "ctb", email = "ebey@uw.edu"),
+             person("Kevin", "Weiss", role = "aut", email = "kevin.weiss@emory.edu"))
 Maintainer: Samuel M. Jenness <samuel.m.jenness@emory.edu>
 Description: EpiModelHIV provides extensions to our general EpiModel package to allow for simulating HIV transmission
   dynamics among two populations: men who have sex with men (MSM) in the United States and heterosexual adults in
@@ -14,20 +15,19 @@ License: GPL-3
 Depends:
     R (>= 3.2.0),
     EpiModel (>= 1.2.7),
+    EpiModelHPC (>= 1.3.1),
     ergm (>= 3.5),
     tergm,
     tergmLite
 Imports:
     bindata,
     network,
     networkDynamic,
-    Rcpp,
     dplyr
 Suggests:
     knitr,
     testthat
 VignetteBuilder: knitr
-LinkingTo: ergm,
-    Rcpp
+LinkingTo: ergm
 RoxygenNote: 6.0.1
 LazyData: true

NAMESPACE

@@ -46,6 +46,9 @@ export(setBirthAttr_het)
 export(simnet_het)
 export(simnet_msm)
 export(sourceDir)
+export(sti_recov)
+export(sti_trans)
+export(sti_tx)
 export(test_msm)
 export(trans_het)
 export(trans_msm)
@@ -59,13 +62,13 @@ export(verbose_msm)
 export(vl_het)
 export(vl_msm)
 import(EpiModel)
+import(EpiModelHPC)
 import(bindata)
 import(ergm)
 import(network)
 import(networkDynamic)
 import(tergm)
 import(tergmLite)
-importFrom(Rcpp,sourceCpp)
 importFrom(dplyr,group_by)
 importFrom(dplyr,summarise)
 importFrom(stats,plogis)

R/EpiModelHIV-package.R

@@ -1,10 +1,10 @@
-#' HIV Transmission Dynamics among MSM and Heterosexuals
+#' Network-Based Epidemic Modeling of HIV Transmission among MSM and Heterosexual Populations
 #'
 #' \tabular{ll}{
 #'    Package: \tab EpiModelHIV\cr
 #'    Type: \tab Package\cr
-#'    Version: \tab 1.0.0\cr
-#'    Date: \tab 2016-06-25\cr
+#'    Version: \tab 1.5.0\cr
+#'    Date: \tab 2017-05-04\cr
 #'    License: \tab GPL-3\cr
 #'    LazyLoad: \tab yes\cr
 #' }
@@ -18,9 +18,8 @@
 #' @name EpiModelHIV-package
 #' @aliases EpiModelHIV
 #'
-#' @import EpiModel network networkDynamic tergmLite tergm ergm bindata
+#' @import EpiModel EpiModelHPC network networkDynamic tergmLite tergm ergm bindata
 #' @importFrom stats rbinom rgeom rmultinom rpois runif simulate rnbinom plogis
-#' @importFrom Rcpp sourceCpp
 #' @importFrom dplyr group_by summarise
 #'
 #' @useDynLib EpiModelHIV

R/mod.acts.R

@@ -82,6 +82,11 @@ acts_msm <- function(dat, at) {
                  (num.B == 0) * base.ai.WW.rate
       ai.rate <- ai.rate * ai.scale
 
+      ## STI associated cessation of activity
+      idsCease <- which(dat$attr$GC.cease == 1 | dat$attr$CT.cease == 1)
+      noActs <- el[, "p1"] %in% idsCease | el[, "p2"] %in% idsCease
+      ai.rate[noActs] <- 0
+
       # Final act number
       if (fixed == FALSE) {
         ai <- rpois(length(ai.rate), ai.rate)

R/RcppExports.R → R/mod.aging.R

@@ -1,9 +1,7 @@
-# Generated by using Rcpp::compileAttributes() -> do not edit by hand
-# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 #' @title Aging Module
 #'
-#' @description Module for aging over time for nodes in the population.
+#' @description Module for aging over time for active nodes in the population.
 #'
 #' @param dat Master data list object of class \code{dat} containing networks,
 #'        individual-level attributes, and summary statistics.
@@ -16,13 +14,26 @@
 #'
 #' @keywords module msm
 #' @export
+#'
 aging_msm <- function(dat, at) {
-    .Call('EpiModelHIV_aging_msm', PACKAGE = 'EpiModelHIV', dat, at)
+
+  time.unit <- dat$param$time.unit
+
+  age <- dat$attr$age
+  active <- dat$attr$active
+
+  age[active == 1] <- age[active == 1] + time.unit / 365
+
+  dat$attr$age <- age
+  dat$attr$sqrt.age <- sqrt(age)
+
+  return(dat)
 }
 
+
 #' @title Aging Module
 #'
-#' @description This module ages all nodes in the population by one time
+#' @description This module ages all active nodes in the population by one time
 #'              unit at each time step.
 #'
 #' @param dat Master data list object of class \code{dat} containing networks,
@@ -31,7 +42,21 @@ aging_msm <- function(dat, at) {
 #'
 #' @keywords module het
 #' @export
+#'
 aging_het <- function(dat, at) {
-    .Call('EpiModelHIV_aging_het', PACKAGE = 'EpiModelHIV', dat, at)
-}
 
+  ## Parameters
+  time.unit <- dat$param$time.unit
+
+  ## Attributes
+  age <- dat$attr$age
+  active <- dat$attr$active
+
+  ## Updates
+  age[active == 1] <- age[active == 1] + time.unit/365
+
+  ## Save out
+  dat$attr$age <- age
+
+  return(dat)
+}

R/mod.condoms.R

@@ -22,15 +22,25 @@
 #'
 condoms_msm <- function(dat, at) {
 
+  # Attributes
+  uid <- dat$attr$uid
+  diag.status <- dat$attr$diag.status
+  race <- dat$attr$race
+  prepStat <- dat$attr$prepStat
+  prepClass <- dat$attr$prepClass
+
+  # Parameters
+  rcomp.prob <- dat$param$rcomp.prob
+  rcomp.adh.groups <- dat$param$rcomp.adh.groups
+  rcomp.main.only <- dat$param$rcomp.main.only
+  rcomp.discl.only <- dat$param$rcomp.discl.only
+
+  el <- dat$temp$el
+
   for (type in c("main", "pers", "inst")) {
 
     ## Variables ##
 
-    # Attributes
-    uid <- dat$attr$uid
-    diag.status <- dat$attr$diag.status
-    race <- dat$attr$race
-
     # Parameters
     cond.rr.BB <- dat$param$cond.rr.BB
     cond.rr.BW <- dat$param$cond.rr.BW
@@ -64,7 +74,6 @@ condoms_msm <- function(dat, at) {
       ptype <- 3
     }
 
-    el <- dat$temp$el
     elt <- el[el[, "ptype"] == ptype, ]
 
     ## Process ##
@@ -77,6 +86,7 @@ condoms_msm <- function(dat, at) {
                  (num.B == 1) * (cond.BW.prob * cond.rr.BW) +
                  (num.B == 0) * (cond.WW.prob * cond.rr.WW)
 
+
     # Transform to UAI logit
     uai.prob <- 1 - cond.prob
     uai.logodds <- log(uai.prob / (1 - uai.prob))
@@ -113,32 +123,56 @@ condoms_msm <- function(dat, at) {
       ca2 <- cond.always[elt[, 2]]
       uai.prob <- ifelse(ca1 == 1 | ca2 == 1, 0, uai.prob)
       if (type == "pers") {
-        dat$epi$cprob.always.pers[at] <- mean(uai.prob == 0)
+        dat$epi$cprob.always.pers <- NULL
+        # dat$epi$cprob.always.pers[at] <- mean(uai.prob == 0)
       } else {
-        dat$epi$cprob.always.inst[at] <- mean(uai.prob == 0)
+        dat$epi$cprob.always.inst <- NULL
+        # dat$epi$cprob.always.inst[at] <- mean(uai.prob == 0)
+      }
+    }
+
+    # PrEP Status (risk compensation)
+    if (rcomp.prob > 0) {
+
+      idsRC <- which((prepStat[elt[, 1]] == 1 & prepClass[elt[, 1]] %in% rcomp.adh.groups) |
+                       (prepStat[elt[, 2]] == 1 & prepClass[elt[, 2]] %in% rcomp.adh.groups))
+
+      if (rcomp.main.only == TRUE & ptype > 1) {
+        idsRC <- NULL
+      }
+      if (rcomp.discl.only == TRUE) {
+        idsRC <- intersect(idsRC, isDisc)
       }
+      uai.prob[idsRC] <- 1 - (1 - uai.prob[idsRC]) * (1 - rcomp.prob)
     }
 
     ai.vec <- elt[, "ai"]
-    pos <- rep(elt[, "p1"], ai.vec)
-    neg <- rep(elt[, "p2"], ai.vec)
+    p1 <- rep(elt[, "p1"], ai.vec)
+    p2 <- rep(elt[, "p2"], ai.vec)
     ptype <- rep(elt[, "ptype"], ai.vec)
 
     uai.prob.peract <- rep(uai.prob, ai.vec)
-    uai <- rbinom(length(pos), 1, uai.prob.peract)
+    uai <- rbinom(length(p1), 1, uai.prob.peract)
 
     if (type == "main") {
       pid <- rep(1:length(ai.vec), ai.vec)
-      al <- cbind(pos, neg, ptype, uai, pid)
+      al <- cbind(p1, p2, ptype, uai, pid)
     } else {
       pid <- rep(max(al[, "pid"]) + (1:length(ai.vec)), ai.vec)
-      tmp.al <- cbind(pos, neg, ptype, uai, pid)
+      tmp.al <- cbind(p1, p2, ptype, uai, pid)
       al <- rbind(al, tmp.al)
     }
 
   } # end ptype loop
 
   dat$temp$al <- al
 
+  if (at == 2) {
+    dat$epi$ai.events <- rep(NA, 2)
+    dat$epi$uai.events <- rep(NA, 2)
+  }
+  dat$epi$ai.events[at] <- nrow(al)
+  dat$epi$uai.events[at] <- sum(al[, "uai"])
+
   return(dat)
 }

R/mod.initialize.R

@@ -41,7 +41,7 @@ initialize_msm <- function(x, param, init, control, s) {
     nw[[i]] <- remove_bad_roles_msm(nw[[i]])
   }
 
-  ## ergm_prep here
+  ## Build initial edgelists
   dat$el <- list()
   dat$p <- list()
   for (i in 1:2) {
@@ -114,6 +114,9 @@ initialize_msm <- function(x, param, init, control, s) {
   dat$attr$prepClass <- rep(NA, num)
   dat$attr$prepElig <- rep(NA, num)
   dat$attr$prepStat <- rep(0, num)
+  dat$attr$prepStartTime <- rep(NA, num)
+  dat$attr$prepLastRisk <- rep(NA, num)
+  dat$attr$prepLastStiScreen <- rep(NA, num)
 
   # One-off AI class
   inst.ai.class <- rep(NA, num)
@@ -136,6 +139,68 @@ initialize_msm <- function(x, param, init, control, s) {
   # HIV-related attributes
   dat <- init_status_msm(dat)
 
+  ## GC/CT status
+  idsUreth <- which(role.class %in% c("I", "V"))
+  idsRect <- which(role.class %in% c("R", "V"))
+
+  uGC <- rGC <- rep(0, num)
+  uCT <- rCT <- rep(0, num)
+
+  # Initialize GC infection at both sites
+  idsUGC <- sample(idsUreth, size = round(init$prev.ugc * num), FALSE)
+  uGC[idsUGC] <- 1
+
+  idsRGC <- sample(setdiff(idsRect, idsUGC), size = round(init$prev.rgc * num), FALSE)
+  rGC[idsRGC] <- 1
+
+  dat$attr$rGC <- rGC
+  dat$attr$uGC <- uGC
+
+  dat$attr$rGC.sympt <- dat$attr$uGC.sympt <- rep(NA, num)
+  dat$attr$rGC.sympt[rGC == 1] <- rbinom(sum(rGC == 1), 1, dat$param$rgc.sympt.prob)
+  dat$attr$uGC.sympt[uGC == 1] <- rbinom(sum(uGC == 1), 1, dat$param$ugc.sympt.prob)
+
+  dat$attr$rGC.infTime <- dat$attr$uGC.infTime <- rep(NA, length(dat$attr$active))
+  dat$attr$rGC.infTime[rGC == 1] <- 1
+  dat$attr$uGC.infTime[uGC == 1] <- 1
+
+  dat$attr$rGC.timesInf <- rep(0, num)
+  dat$attr$rGC.timesInf[rGC == 1] <- 1
+  dat$attr$uGC.timesInf <- rep(0, num)
+  dat$attr$uGC.timesInf[uGC == 1] <- 1
+
+  dat$attr$rGC.tx <- dat$attr$uGC.tx <- rep(NA, num)
+  dat$attr$rGC.tx.prep <- dat$attr$uGC.tx.prep <- rep(NA, num)
+  dat$attr$GC.cease <- rep(NA, num)
+
+  # Initialize CT infection at both sites
+  idsUCT <- sample(idsUreth, size = round(init$prev.uct * num), FALSE)
+  uCT[idsUCT] <- 1
+
+  idsRCT <- sample(setdiff(idsRect, idsUCT), size = round(init$prev.rct * num), FALSE)
+  rCT[idsRCT] <- 1
+
+  dat$attr$rCT <- rCT
+  dat$attr$uCT <- uCT
+
+  dat$attr$rCT.sympt <- dat$attr$uCT.sympt <- rep(NA, num)
+  dat$attr$rCT.sympt[rCT == 1] <- rbinom(sum(rCT == 1), 1, dat$param$rct.sympt.prob)
+  dat$attr$uCT.sympt[uCT == 1] <- rbinom(sum(uCT == 1), 1, dat$param$uct.sympt.prob)
+
+  dat$attr$rCT.infTime <- dat$attr$uCT.infTime <- rep(NA, num)
+  dat$attr$rCT.infTime[dat$attr$rCT == 1] <- 1
+  dat$attr$uCT.infTime[dat$attr$uCT == 1] <- 1
+
+  dat$attr$rCT.timesInf <- rep(0, num)
+  dat$attr$rCT.timesInf[rCT == 1] <- 1
+  dat$attr$uCT.timesInf <- rep(0, num)
+  dat$attr$uCT.timesInf[uCT == 1] <- 1
+
+  dat$attr$rCT.tx <- dat$attr$uCT.tx <- rep(NA, num)
+  dat$attr$rCT.tx.prep <- dat$attr$uCT.tx.prep <- rep(NA, num)
+  dat$attr$CT.cease <- rep(NA, num)
+
+
   # CCR5
   dat <- init_ccr5_msm(dat)
 
@@ -149,12 +214,6 @@ initialize_msm <- function(x, param, init, control, s) {
   dat$temp$discl.list <- matrix(NA, nrow = 0, ncol = 3)
   colnames(dat$temp$discl.list) <- c("pos", "neg", "discl.time")
 
-  if (control$save.nwstats == TRUE) {
-    dat$stats <- list()
-    dat$stats$nwstats <- list()
-
-  }
-
   dat <- prevalence_msm(dat, at = 1)
 
   class(dat) <- "dat"