From e319c7951a72feaa57ff64da8c375ef9a25e9802 Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Mon, 8 Aug 2016 23:40:55 -0400
Subject: [PATCH 01/40] Port over STI model

---
 NAMESPACE          |   3 +
 R/mod.acts.R       |   5 +
 R/mod.condoms.R    |  60 ++++-
 R/mod.initialize.R |  70 +++++-
 R/mod.position.R   |  41 ++--
 R/mod.prep.R       |  88 +++----
 R/mod.prevalence.R | 100 ++++++--
 R/mod.riskhist.R   |  75 ++----
 R/mod.sti.R        | 589 +++++++++++++++++++++++++++++++++++++++++++++
 R/mod.trans.R      |  84 +++++--
 R/mod.verbose.R    |  27 ++-
 R/params.R         | 114 ++++++++-
 man/control_msm.Rd |  13 +-
 man/param_msm.Rd   | 105 +++++++-
 man/sti_recov.Rd   |  20 ++
 man/sti_trans.Rd   |  21 ++
 man/sti_tx.Rd      |  20 ++
 17 files changed, 1239 insertions(+), 196 deletions(-)
 create mode 100644 R/mod.sti.R
 create mode 100644 man/sti_recov.Rd
 create mode 100644 man/sti_trans.Rd
 create mode 100644 man/sti_tx.Rd

diff --git a/NAMESPACE b/NAMESPACE
index b706ae2b..d8e0a138 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/mod.acts.R b/R/mod.acts.R
index 2ed0f42c..4ab27b41 100644
--- a/R/mod.acts.R
+++ b/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)
diff --git a/R/mod.condoms.R b/R/mod.condoms.R
index 6cd8607a..0a123b42 100644
--- a/R/mod.condoms.R
+++ b/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,26 +123,43 @@ 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)
     }
 
@@ -140,5 +167,12 @@ condoms_msm <- function(dat, at) {
 
   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)
 }
diff --git a/R/mod.initialize.R b/R/mod.initialize.R
index f45ab29f..ce974408 100644
--- a/R/mod.initialize.R
+++ b/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,65 @@ 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$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$CT.cease <- rep(NA, num)
+
   # CCR5
   dat <- init_ccr5_msm(dat)
 
@@ -149,12 +211,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"
diff --git a/R/mod.position.R b/R/mod.position.R
index 0dc48dc7..5c24aea5 100644
--- a/R/mod.position.R
+++ b/R/mod.position.R
@@ -20,7 +20,7 @@
 #' susceptible node is insertive for that act.
 #'
 #' @keywords module msm
-#' 
+#'
 #' @export
 #'
 position_msm <- function(dat, at) {
@@ -31,49 +31,50 @@ position_msm <- function(dat, at) {
     return(dat)
   }
 
-  status <- dat$attr$status
-  dal <- al[which(status[al[, 1]] == 1 & status[al[, 2]] == 0), ]
-  dat$temp$al <- NULL
+  # Attributes
 
   role.class <- dat$attr$role.class
   ins.quot <- dat$attr$ins.quot
   race <- dat$attr$race
 
+  # Parameters
+
   vv.iev.BB.prob <- dat$param$vv.iev.BB.prob
   vv.iev.BW.prob <- dat$param$vv.iev.BW.prob
   vv.iev.WW.prob <- dat$param$vv.iev.WW.prob
 
 
   ## Process
-  pos.role.class <- role.class[dal[, 1]]
-  neg.role.class <- role.class[dal[, 2]]
+  p1.role.class <- role.class[al[, 1]]
+  p2.role.class <- role.class[al[, 2]]
 
-  ins <- rep(NA, length(pos.role.class))
-  ins[which(pos.role.class == "I")] <- 1  # "P"
-  ins[which(pos.role.class == "R")] <- 0  # "N"
-  ins[which(neg.role.class == "I")] <- 0  # "N"
-  ins[which(neg.role.class == "R")] <- 1  # "P"
+  ins <- rep(NA, length(p1.role.class))
+  ins[which(p1.role.class == "I")] <- 1
+  ins[which(p1.role.class == "R")] <- 0
+  ins[which(p2.role.class == "I")] <- 0
+  ins[which(p2.role.class == "R")] <- 1
 
-  vv <- which(pos.role.class == "V" & neg.role.class == "V")
-  vv.race.combo <- paste0(race[dal[, 1]][vv], race[dal[, 2]][vv])
+  vv <- which(p1.role.class == "V" & p2.role.class == "V")
+  vv.race.combo <- paste0(race[al[, 1]][vv], race[al[, 2]][vv])
   vv.race.combo[vv.race.combo == "WB"] <- "BW"
   vv.iev.prob <- (vv.race.combo == "BB") * vv.iev.BB.prob +
                  (vv.race.combo == "BW") * vv.iev.BW.prob +
                  (vv.race.combo == "WW") * vv.iev.WW.prob
 
+  # intra-event versatility
   iev <- rbinom(length(vv), 1, vv.iev.prob)
-  ins[vv[iev == 1]] <- 2 # "B"
+  ins[vv[iev == 1]] <- 2 # both are insertive, acts will be doubled
   vv.remaining <- vv[iev == 0]
 
-  inspos.prob <- ins.quot[dal[, 1][vv.remaining]] /
-                 (ins.quot[dal[, 1][vv.remaining]] + ins.quot[dal[, 2][vv.remaining]])
-  inspos <- rbinom(length(vv.remaining), 1, inspos.prob)
-  ins[vv.remaining[inspos == 1]] <- 1  # "P"
-  ins[vv.remaining[inspos == 0]] <- 0  # "N"
+  p1.ins.prob <- ins.quot[al[, 1][vv.remaining]] /
+    (ins.quot[al[, 1][vv.remaining]] + ins.quot[al[, 2][vv.remaining]])
+  p1.ins <- rbinom(length(vv.remaining), 1, p1.ins.prob)
+  ins[vv.remaining[p1.ins == 1]] <- 1
+  ins[vv.remaining[p1.ins == 0]] <- 0
 
 
   ## Output
-  dat$temp$dal <- cbind(dal, ins)
+  dat$temp$al <- cbind(al, ins)
 
   return(dat)
 }
diff --git a/R/mod.prep.R b/R/mod.prep.R
index 41e5ce5c..44afed97 100644
--- a/R/mod.prep.R
+++ b/R/mod.prep.R
@@ -17,83 +17,67 @@ prep_msm <- function(dat, at) {
   }
 
   ## Variables
+
+  # Attributes
+  active <- dat$attr$active
   status <- dat$attr$status
   diag.status <- dat$attr$diag.status
   lnt <- dat$attr$last.neg.test
-
   prepElig <- dat$attr$prepElig
   prepStat <- dat$attr$prepStat
   prepClass <- dat$attr$prepClass
+  prepLastRisk <- dat$attr$prepLastRisk
+  prepStartTime <- dat$attr$prepStartTime
+
+  # Parameters
 
-  prep.elig.model <- dat$param$prep.elig.model
   prep.coverage <- dat$param$prep.coverage
   prep.cov.rate <- dat$param$prep.cov.rate
   prep.class.prob <- dat$param$prep.class.prob
-  prep.risk.reassess <- dat$param$prep.risk.reassess
 
 
   ## Eligibility ---------------------------------------------------------------
 
   # Base eligibility
-  idsEligStart <- which(status == 0 & prepStat == 0 & lnt == at)
+  idsEligStart <- which(active == 1 & status == 0 & prepStat == 0 & lnt == at)
 
-  idsEligStop <- NULL
-  if (prep.risk.reassess == TRUE) {
-    idsEligStop <- which(prepStat == 1 & lnt == at)
-  }
+  # Core eligiblity
+  ind1 <- dat$attr$prep.ind.uai.mono
+  ind2 <- dat$attr$prep.ind.uai.nmain
+  ind3 <- dat$attr$prep.ind.ai.sd
+  ind4 <- dat$attr$prep.ind.sti
 
   twind <- at - dat$param$prep.risk.int
-
-  # Core eligiblity scenarios
-  if (prep.elig.model != "base") {
-    if (substr(prep.elig.model, 1, 3) == "cdc") {
-      if (prep.elig.model == "cdc1") {
-        c1 <- dat$attr$uai.mono2.nt.6mo
-        c2 <- dat$attr$uai.nonmonog
-        c3 <- dat$attr$ai.sd.mc
-      } else if (prep.elig.model == "cdc2") {
-        c1 <- dat$attr$uai.mono2.nt.6mo
-        c2 <- dat$attr$uai.nmain
-        c3 <- dat$attr$ai.sd.mc
-      } else if (prep.elig.model == "cdc3") {
-        c1 <- dat$attr$uai.mono1.nt.6mo
-        c2 <- dat$attr$uai.nmain
-        c3 <- dat$attr$ai.sd.mc
-      } else if (prep.elig.model == "cdc4") {
-        c1 <- dat$attr$uai.mono1.nt.6mo
-        c2 <- dat$attr$uai.nmain
-        c3 <- dat$attr$uai.sd.mc
-      }
-      idsEligStart <- intersect(which(c1 >= twind |
-                                      c2 >= twind |
-                                      c3 >= twind),
-                                idsEligStart)
-      idsEligStop <- intersect(which(c1 < twind &
-                                     c2 < twind &
-                                     c3 < twind),
-                                idsEligStop)
-    } else {
-      c1 <- dat$attr[[prep.elig.model]]
-      idsEligStart <- intersect(which(c1 >= twind), idsEligStart)
-      idsEligStop <- intersect(which(c1 < twind), idsEligStop)
-    }
-  }
+  idsEligStart <- intersect(which(ind1 >= twind | ind2 >= twind |
+                                    ind3 >= twind | ind4 >= twind),
+                            idsEligStart)
 
   prepElig[idsEligStart] <- 1
-  prepElig[idsEligStop] <- 0
 
 
   ## Stoppage ------------------------------------------------------------------
 
+  # No indications
+  idsRiskAssess <- which(active == 1 & prepStat == 1 & lnt == at & (at - prepLastRisk) >= 52)
+  prepLastRisk[idsRiskAssess] <- at
+
+  idsEligStop <- intersect(which(ind1 < twind & ind2 < twind &
+                                   ind3 < twind & ind4 < twind),
+                           idsRiskAssess)
+
+  prepElig[idsEligStop] <- 0
+
   # Diagnosis
-  idsStpDx <- which(prepStat == 1 & diag.status == 1)
+  idsStpDx <- which(active == 1 & prepStat == 1 & diag.status == 1)
 
-  # Transition to ineligibility
-  idsStpInelig <- idsEligStop
+  # Death
+  idsStpDth <- which(active == 0 & prepStat == 1)
 
   # Reset PrEP status
-  idsStp <- c(idsStpDx, idsStpInelig)
+  idsStp <- c(idsStpDx, idsStpDth, idsEligStop)
   prepStat[idsStp] <- 0
+  prepLastRisk[idsStp] <- NA
+  prepStartTime[idsStp] <- NA
 
 
   ## Initiation ----------------------------------------------------------------
@@ -105,7 +89,7 @@ prep_msm <- function(dat, at) {
   nEligSt <- length(idsEligSt)
 
   nStart <- max(0, min(nEligSt, round((prep.coverage - prepCov) *
-                                             sum(prepElig == 1, na.rm = TRUE))))
+                                        sum(prepElig == 1, na.rm = TRUE))))
   idsStart <- NULL
   if (nStart > 0) {
     if (prep.cov.rate >= 1) {
@@ -118,8 +102,10 @@ prep_msm <- function(dat, at) {
   # Attributes
   if (length(idsStart) > 0) {
     prepStat[idsStart] <- 1
+    prepStartTime[idsStart] <- at
+    prepLastRisk[idsStart] <- at
 
-    # PrEP class is fixed over PrEP cycles
+    # PrEP class
     needPC <- which(is.na(prepClass[idsStart]))
     prepClass[idsStart[needPC]] <- sample(x = 0:3, size = length(needPC),
                                           replace = TRUE, prob = prep.class.prob)
@@ -131,7 +117,9 @@ prep_msm <- function(dat, at) {
   # Attributes
   dat$attr$prepElig <- prepElig
   dat$attr$prepStat <- prepStat
+  dat$attr$prepStartTime <- prepStartTime
   dat$attr$prepClass <- prepClass
+  dat$attr$prepLastRisk <- prepLastRisk
 
   # Summary Statistics
   dat$epi$prepCov[at] <- prepCov
diff --git a/R/mod.prevalence.R b/R/mod.prevalence.R
index f2e4d00c..3ea3e8fe 100644
--- a/R/mod.prevalence.R
+++ b/R/mod.prevalence.R
@@ -24,9 +24,24 @@
 #'
 prevalence_msm <- function(dat, at) {
 
+  ## Variables
+
+  # Attributes
+
+  active <- dat$attr$active
   race <- dat$attr$race
   status <- dat$attr$status
   prepStat <- dat$attr$prepStat
+  prepElig <- dat$attr$prepElig
+  rGC <- dat$attr$rGC
+  uGC <- dat$attr$uGC
+  rCT <- dat$attr$rCT
+  uCT <- dat$attr$uCT
+  rGC.sympt <- dat$attr$rGC.sympt
+  uGC.sympt <- dat$attr$uGC.sympt
+  rCT.sympt <- dat$attr$rCT.sympt
+  uCT.sympt <- dat$attr$uCT.sympt
+
 
   nsteps <- dat$control$nsteps
   rNA <- rep(NA, nsteps)
@@ -42,26 +57,59 @@ prevalence_msm <- function(dat, at) {
     dat$epi$i.prev <- rNA
     dat$epi$i.prev.B <- rNA
     dat$epi$i.prev.W <- rNA
-    dat$epi$nBirths <- rNA
-    dat$epi$dth.gen <- rNA
-    dat$epi$dth.dis <- rNA
     dat$epi$incid <- rNA
+    dat$epi$ir100 <- rNA
 
     dat$epi$prepCurr <- rNA
     dat$epi$prepCov <- rNA
     dat$epi$prepElig <- rNA
     dat$epi$prepStart <- rNA
-    dat$epi$incid.prep0 <- rNA
-    dat$epi$incid.prep1 <- rNA
     dat$epi$i.num.prep0 <- rNA
     dat$epi$i.num.prep1 <- rNA
 
-    dat$epi$cprob.always.pers <- rNA
-    dat$epi$cprob.always.inst <- rNA
+    dat$epi$prev.rgc <- rNA
+    dat$epi$prev.ugc <- rNA
+    dat$epi$prev.gc <- rNA
+    dat$epi$prev.gc.sympt <- rNA
+    dat$epi$prev.gc.dual <- rNA
+
+    dat$epi$prev.rct <- rNA
+    dat$epi$prev.uct <- rNA
+    dat$epi$prev.ct <- rNA
+    dat$epi$prev.ct.sympt <- rNA
+    dat$epi$prev.ct.dual <- rNA
+
+    dat$epi$prev.rgcct <- rNA
+    dat$epi$prev.ugcct <- rNA
+
+    dat$epi$incid.rgc <- rNA
+    dat$epi$incid.ugc <- rNA
+    dat$epi$incid.gc <- rNA
+    dat$epi$incid.rct <- rNA
+    dat$epi$incid.uct <- rNA
+    dat$epi$incid.ct <- rNA
+
+    dat$epi$ir100.rgc <- rNA
+    dat$epi$ir100.ugc <- rNA
+    dat$epi$ir100.gc <- rNA
+    dat$epi$ir100.rct <- rNA
+    dat$epi$ir100.uct <- rNA
+    dat$epi$ir100.ct <- rNA
+
+    dat$epi$recov.rgc <- rNA
+    dat$epi$recov.ugc <- rNA
+    dat$epi$recov.rct <- rNA
+    dat$epi$recov.uct <- rNA
+
+    dat$epi$trans.main <- rNA
+    dat$epi$trans.casl <- rNA
+    dat$epi$trans.inst <- rNA
+
+    dat$epi$txGC <- rNA
+    dat$epi$txCT <- rNA
   }
 
-
-  dat$epi$num[at] <- length(status)
+  dat$epi$num[at] <- sum(active == 1, na.rm = TRUE)
   dat$epi$num.B[at] <- sum(race == "B", na.rm = TRUE)
   dat$epi$num.W[at] <- sum(race == "W", na.rm = TRUE)
   dat$epi$s.num[at] <- sum(status == 0, na.rm = TRUE)
@@ -71,21 +119,43 @@ prevalence_msm <- function(dat, at) {
   dat$epi$i.prev[at] <- dat$epi$i.num[at] / dat$epi$num[at]
   dat$epi$i.prev.B[at] <- dat$epi$i.num.B[at] / dat$epi$num.B[at]
   dat$epi$i.prev.W[at] <- dat$epi$i.num.W[at] / dat$epi$num.W[at]
+  dat$epi$ir100[at] <- (dat$epi$incid[at] / sum(status == 0, na.rm = TRUE)) * 5200
 
   dat$epi$prepCurr[at] <- sum(prepStat == 1, na.rm = TRUE)
-  dat$epi$prepElig[at] <- sum(dat$attr$prepElig == 1, na.rm = TRUE)
-  dat$epi$i.num.prep0[at] <- sum((is.na(prepStat) | prepStat == 0) &
-                                 status == 1, na.rm = TRUE)
+  dat$epi$prepElig[at] <- sum(prepElig == 1, na.rm = TRUE)
+  dat$epi$i.num.prep0[at] <- sum((is.na(prepStat) | prepStat == 0) & status == 1, na.rm = TRUE)
   dat$epi$i.num.prep1[at] <- sum(prepStat == 1 & status == 1, na.rm = TRUE)
   dat$epi$i.prev.prep0[at] <- dat$epi$i.num.prep0[at] /
-                              sum((is.na(prepStat) | prepStat == 0), na.rm = TRUE)
+    sum((is.na(prepStat) | prepStat == 0), na.rm = TRUE)
   if (at == 1) {
     dat$epi$i.prev.prep1[1] <- 0
   } else {
-    dat$epi$i.prev.prep1[at] <- dat$epi$i.num.prep1[at] /
-                                sum(prepStat == 1, na.rm = TRUE)
+    dat$epi$i.prev.prep1[at] <- dat$epi$i.num.prep1[at] / sum(prepStat == 1, na.rm = TRUE)
   }
 
+  dat$epi$prev.rgc[at] <- sum(rGC == 1, na.rm = TRUE) / dat$epi$num[at]
+  dat$epi$prev.ugc[at] <- sum(uGC == 1, na.rm = TRUE) / dat$epi$num[at]
+  dat$epi$prev.gc[at] <- sum((rGC == 1 | uGC == 1), na.rm = TRUE) / dat$epi$num[at]
+  dat$epi$prev.gc.sympt[at] <- sum((rGC.sympt == 1 | uGC.sympt == 1)) / dat$epi$num[at]
+  dat$epi$prev.gc.dual[at] <- sum((rGC == 1 & uGC == 1), na.rm = TRUE) / dat$epi$num[at]
+
+  dat$epi$prev.rct[at] <- sum(rCT == 1, na.rm = TRUE) / dat$epi$num[at]
+  dat$epi$prev.uct[at] <- sum(uCT == 1, na.rm = TRUE) / dat$epi$num[at]
+  dat$epi$prev.ct[at] <- sum((rCT == 1 | uCT == 1), na.rm = TRUE) / dat$epi$num[at]
+  dat$epi$prev.ct.sympt[at] <- sum((rCT.sympt == 1 | uCT.sympt == 1)) / dat$epi$num[at]
+  dat$epi$prev.ct.dual[at] <- sum((rCT == 1 & uCT == 1), na.rm = TRUE) / dat$epi$num[at]
+
+  dat$epi$prev.rgcct[at] <- sum(rGC == 1 | rCT == 1, na.rm = TRUE) / dat$epi$num[at]
+  dat$epi$prev.ugcct[at] <- sum(uGC == 1 | uCT == 1, na.rm = TRUE) / dat$epi$num[at]
+
+  dat$epi$ir100.rgc[at] <- (dat$epi$incid.rgc[at] / sum(rGC == 0, na.rm = TRUE)) * 5200
+  dat$epi$ir100.ugc[at] <- (dat$epi$incid.ugc[at] / sum(uGC == 0, na.rm = TRUE)) * 5200
+  dat$epi$ir100.gc[at] <- (dat$epi$incid.gc[at]/ sum(rGC == 0 | uGC == 0, na.rm = TRUE)) * 5200
+
+  dat$epi$ir100.rct[at] <- (dat$epi$incid.rct[at] / sum(rCT == 0, na.rm = TRUE)) * 5200
+  dat$epi$ir100.uct[at] <- (dat$epi$incid.uct[at] / sum(uCT == 0, na.rm = TRUE)) * 5200
+  dat$epi$ir100.ct[at] <- (dat$epi$incid.ct[at]/ sum(rCT == 0 | uCT == 0, na.rm = TRUE)) * 5200
+
   return(dat)
 }
 
diff --git a/R/mod.riskhist.R b/R/mod.riskhist.R
index dd91fbab..e6978070 100644
--- a/R/mod.riskhist.R
+++ b/R/mod.riskhist.R
@@ -18,8 +18,15 @@ riskhist_msm <- function(dat, at) {
 
   ## Attributes
   uid <- dat$attr$uid
+  dx <- dat$attr$diag.status
+  since.test <- at - dat$attr$last.neg.test
+  rGC.tx <- dat$attr$rGC.tx
+  uGC.tx <- dat$attr$uGC.tx
+  rCT.tx <- dat$attr$rCT.tx
+  uCT.tx <- dat$attr$uCT.tx
 
   ## Parameters
+  time.unit <- dat$param$time.unit
 
   ## Edgelist, adds uai summation per partnership from act list
   pid <- NULL # For R CMD Check
@@ -31,13 +38,12 @@ riskhist_msm <- function(dat, at) {
   # Remove concordant positive edges
   el2 <- el[el$st2 == 0, ]
 
-  if (is.null(dat$attr$uai.mono2.nt.3mo)) {
-    dat$attr$uai.mono2.nt.3mo <- rep(NA, length(uid))
-    dat$attr$uai.mono1.nt.3mo <- rep(NA, length(uid))
-    dat$attr$uai.nonmonog <- rep(NA, length(uid))
-    dat$attr$uai.nmain <- rep(NA, length(uid))
-    dat$attr$ai.sd.mc <- rep(NA, length(uid))
-    dat$attr$uai.sd.mc <- rep(NA, length(uid))
+  # Initialize attributes
+  if (is.null(dat$attr$prep.ind.uai.mono)) {
+    dat$attr$prep.ind.uai.mono <- rep(NA, length(uid))
+    dat$attr$prep.ind.uai.nmain <- rep(NA, length(uid))
+    dat$attr$prep.ind.ai.sd <- rep(NA, length(uid))
+    dat$attr$prep.ind.sti <- rep(NA, length(uid))
   }
 
   ## Degree ##
@@ -56,57 +62,23 @@ riskhist_msm <- function(dat, at) {
   tot.deg <- main.deg + casl.deg + inst.deg
   uai.mono1 <- intersect(which(tot.deg == 1), uai.any)
 
-  # Monogamous partnerships: 2-sided
-  mono.2s <- tot.deg[el2$p1] == 1 & tot.deg[el2$p2] == 1
-  ai.mono2 <- sort(unname(do.call("c", c(el2[mono.2s, 1:2]))))
-  uai.mono2 <- intersect(ai.mono2, uai.any)
-
   # "Negative" partnerships
   tneg <- unique(c(el2$p1[el2$st1 == 0], el2$p2[el2$st1 == 0]))
-  dx <- dat$attr$diag.status
   fneg <- unique(c(el2$p1[which(dx[el2$p1] == 0)], el2$p2[which(dx[el2$p1] == 0)]))
   all.neg <- c(tneg, fneg)
-  since.test <- at - dat$attr$last.neg.test
-
-
-  ## Condition 1a: UAI in 2-sided monogamous "negative" partnership,
-  ##               partner not tested in past 3, 6 months
-  uai.mono2.neg <- intersect(uai.mono2, all.neg)
-  part.id2 <- c(el2[el2$p1 %in% uai.mono2.neg, 2], el2[el2$p2 %in% uai.mono2.neg, 1])
-  not.tested.3mo <- since.test[part.id2] > (90/dat$param$time.unit)
-  part.not.tested.3mo <- uai.mono2.neg[which(not.tested.3mo == TRUE)]
-  dat$attr$uai.mono2.nt.3mo[part.not.tested.3mo] <- at
-
-  not.tested.6mo <- since.test[part.id2] > (180/dat$param$time.unit)
-  part.not.tested.6mo <- uai.mono2.neg[which(not.tested.6mo == TRUE)]
-  dat$attr$uai.mono2.nt.6mo[part.not.tested.6mo] <- at
-
 
   ## Condition 1b: UAI in 1-sided "monogamous" "negative" partnership,
-  ##               partner not tested in past 3, 6 months
+  ##               partner not tested in past 6 months
   uai.mono1.neg <- intersect(uai.mono1, all.neg)
   part.id1 <- c(el2[el2$p1 %in% uai.mono1.neg, 2], el2[el2$p2 %in% uai.mono1.neg, 1])
-  not.tested.3mo <- since.test[part.id1] > (90/dat$param$time.unit)
-  part.not.tested.3mo <- uai.mono1.neg[which(not.tested.3mo == TRUE)]
-  dat$attr$uai.mono1.nt.3mo[part.not.tested.3mo] <- at
-
-  not.tested.6mo <- since.test[part.id1] > (180/dat$param$time.unit)
+  not.tested.6mo <- since.test[part.id1] > (180/time.unit)
   part.not.tested.6mo <- uai.mono1.neg[which(not.tested.6mo == TRUE)]
-  dat$attr$uai.mono1.nt.6mo[part.not.tested.6mo] <- at
-
-
-  ## Condition 2a: UAI in non-monogamous partnerships
-  el2.uai <- el2[el2$uai > 0, ]
-  vec <- c(el2.uai[, 1], el2.uai[, 2])
-  uai.nonmonog <- unique(vec[duplicated(vec)])
-  dat$attr$uai.nonmonog[uai.nonmonog] <- at
-
+  dat$attr$prep.ind.uai.mono[part.not.tested.6mo] <- at
 
   ## Condition 2b: UAI in non-main partnerships
   uai.nmain <- unique(c(el2$p1[el2$st1 == 0 & el2$uai > 0 & el2$ptype %in% 2:3],
                         el2$p2[el2$uai > 0 & el2$ptype %in% 2:3]))
-  dat$attr$uai.nmain[uai.nmain] <- at
-
+  dat$attr$prep.ind.uai.nmain[uai.nmain] <- at
 
   ## Condition 3a: AI within known serodiscordant partnerships
   el2.cond3 <- el2[el2$st1 == 1 & el2$ptype %in% 1:2, ]
@@ -116,14 +88,13 @@ riskhist_msm <- function(dat, at) {
   disclose.cdl <- discl.list[, 1] * 1e7 + discl.list[, 2]
   delt.cdl <- uid[el2.cond3[, 1]] * 1e7 + uid[el2.cond3[, 2]]
   discl <- (delt.cdl %in% disclose.cdl)
+  ai.sd <- el2.cond3$p2[discl == TRUE]
+  dat$attr$prep.ind.ai.sd[ai.sd] <- at
 
-  ai.sd.mc <- el2.cond3$p2[discl == TRUE]
-  dat$attr$ai.sd.mc[ai.sd.mc] <- at
-
-
-  ## Condition 3b: UAI within known serodiscordant partnerships
-  uai.sd.mc <- el2.cond3$p2[discl == TRUE & el2.cond3$uai > 0]
-  dat$attr$uai.sd.mc[uai.sd.mc] <- at
+  ## Condition 4, any STI diagnosis
+  idsDx <- which(rGC.tx == 1 | uGC.tx == 1 |
+                   rCT.tx == 1 | uCT.tx == 1)
+  dat$attr$prep.ind.sti[idsDx] <- at
 
   return(dat)
 }
diff --git a/R/mod.sti.R b/R/mod.sti.R
new file mode 100644
index 00000000..25fbcf9e
--- /dev/null
+++ b/R/mod.sti.R
@@ -0,0 +1,589 @@
+
+#' @title STI Transmission Module
+#'
+#' @description Stochastically simulates GC/CT transmission given the current
+#'              state of the edgelist.
+#'
+#' @inheritParams aging_msm
+#'
+#' @keywords module msm
+#'
+#' @export
+#'
+sti_trans <- function(dat, at) {
+
+  # Parameters ----------------------------------------------------------
+
+  # Acquisition probabilities given contact with infected man
+  rgc.tprob <- dat$param$rgc.tprob
+  ugc.tprob <- dat$param$ugc.tprob
+  rct.tprob <- dat$param$rct.tprob
+  uct.tprob <- dat$param$uct.tprob
+
+  # Probability of symptoms given infection
+  rgc.sympt.prob <- dat$param$rgc.sympt.prob
+  ugc.sympt.prob <- dat$param$ugc.sympt.prob
+  rct.sympt.prob <- dat$param$rct.sympt.prob
+  uct.sympt.prob <- dat$param$uct.sympt.prob
+
+  # Relative risk of infection given condom use during act
+  sti.cond.rr <- dat$param$sti.cond.rr
+
+  # Cessation
+  gc.prob.cease <- dat$param$gc.prob.cease
+  ct.prob.cease <- dat$param$ct.prob.cease
+
+  # Attributes ----------------------------------------------------------
+
+  # Current infection state
+  rGC <- dat$attr$rGC
+  uGC <- dat$attr$uGC
+  rCT <- dat$attr$rCT
+  uCT <- dat$attr$uCT
+
+  # n Times infected
+  rGC.timesInf <- dat$attr$rGC.timesInf
+  uGC.timesInf <- dat$attr$uGC.timesInf
+  rCT.timesInf <- dat$attr$rCT.timesInf
+  uCT.timesInf <- dat$attr$uCT.timesInf
+
+  # Set disease status to 0 for new births
+  newBirths <- which(dat$attr$arrival.time == at)
+  rGC[newBirths] <- rGC.timesInf[newBirths] <- 0
+  uGC[newBirths] <- uGC.timesInf[newBirths] <- 0
+  rCT[newBirths] <- rCT.timesInf[newBirths] <- 0
+  uCT[newBirths] <- uCT.timesInf[newBirths] <- 0
+
+  # Infection time
+  rGC.infTime <- dat$attr$rGC.infTime
+  uGC.infTime <- dat$attr$uGC.infTime
+  rCT.infTime <- dat$attr$rCT.infTime
+  uCT.infTime <- dat$attr$uCT.infTime
+
+
+
+  # Infection symptoms (non-varying)
+  rGC.sympt <- dat$attr$rGC.sympt
+  uGC.sympt <- dat$attr$uGC.sympt
+  rCT.sympt <- dat$attr$rCT.sympt
+  uCT.sympt <- dat$attr$uCT.sympt
+
+  # Men who cease sexual activity during symptomatic infection
+  GC.cease <- dat$attr$GC.cease
+  CT.cease <- dat$attr$CT.cease
+
+  # Pull act list
+  al <- dat$temp$al
+
+  ## ins variable coding
+  # ins = 0 : p2 is insertive
+  # ins = 1 : p1 is insertive
+  # ins = 2 : both p1 and p2 are insertive
+
+  # Rectal GC -----------------------------------------------------------
+
+  # Requires: uGC in insertive man, and no rGC in receptive man
+  p1Inf_rgc <- which(uGC[al[, "p1"]] == 1 & uGC.infTime[al[, "p1"]] < at &
+                       rGC[al[, "p2"]] == 0 & al[, "ins"] %in% c(1, 2))
+  p2Inf_rgc <- which(uGC[al[, "p1"]] == 0 & uGC.infTime[al[, "p2"]] < at &
+                       rGC[al[, "p1"]] == 0 & al[, "ins"] %in% c(0, 2))
+  allActs_rgc <- c(p1Inf_rgc, p2Inf_rgc)
+
+  # UAI modifier
+  uai_rgc <- al[, "uai"][allActs_rgc]
+  tprob_rgc <- rep(rgc.tprob, length(allActs_rgc))
+  tprob_rgc[uai_rgc == 0] <- tprob_rgc[uai_rgc == 0] * sti.cond.rr
+
+  # Stochastic transmission
+  trans_rgc <- rbinom(length(allActs_rgc), 1, tprob_rgc)
+
+  # Determine the infected partner
+  idsInf_rgc <- NULL
+  if (sum(trans_rgc) > 0) {
+    transAL_rgc <- al[allActs_rgc[trans_rgc == 1], , drop = FALSE]
+    idsInf_rgc <- unique(ifelse(uGC[transAL_rgc[, "p1"]] == 1,
+                                transAL_rgc[, "p2"], transAL_rgc[, "p1"]))
+  }
+
+  # Update attributes
+  rGC[idsInf_rgc] <- 1
+  rGC.infTime[idsInf_rgc] <- at
+  rGC.sympt[idsInf_rgc] <- rbinom(length(idsInf_rgc), 1, rgc.sympt.prob)
+  rGC.timesInf[idsInf_rgc] <- rGC.timesInf[idsInf_rgc] + 1
+
+  # Urethral GC ---------------------------------------------------------
+
+  # Requires: rGC in receptive man, and no uGC in insertive man
+  p1Inf_ugc <- which(rGC[al[, "p1"]] == 1 & rGC.infTime[al[, "p1"]] < at &
+                       uGC[al[, "p2"]] == 0 & al[, "ins"] %in% c(0, 2))
+  p2Inf_ugc <- which(rGC[al[, "p1"]] == 0 & rGC.infTime[al[, "p2"]] < at &
+                       uGC[al[, "p1"]] == 0 & al[, "ins"] %in% c(1, 2))
+  allActs_ugc <- c(p1Inf_ugc, p2Inf_ugc)
+
+  # UAI modifier
+  uai_ugc <- al[, "uai"][allActs_ugc]
+  tprob_ugc <- rep(ugc.tprob, length(allActs_ugc))
+  tprob_ugc[uai_ugc == 0] <- tprob_ugc[uai_ugc == 0] * sti.cond.rr
+
+  # Stochastic transmission
+  trans_ugc <- rbinom(length(allActs_ugc), 1, tprob_ugc)
+
+  # Determine the newly infected partner
+  idsInf_ugc <- NULL
+  if (sum(trans_ugc) > 0) {
+    transAL_ugc <- al[allActs_ugc[trans_ugc == 1],  , drop = FALSE]
+    idsInf_ugc <- unique(ifelse(uGC[transAL_ugc[, "p1"]] == 1,
+                                transAL_ugc[, "p2"], transAL_ugc[, "p1"]))
+  }
+
+  # Update attributes
+  uGC[idsInf_ugc] <- 1
+  uGC.infTime[idsInf_ugc] <- at
+  uGC.sympt[idsInf_ugc] <- rbinom(length(idsInf_ugc), 1, ugc.sympt.prob)
+  uGC.timesInf[idsInf_ugc] <- uGC.timesInf[idsInf_ugc] + 1
+
+
+  # Rectal CT -----------------------------------------------------------
+
+  # Requires: uCT in insertive man, and no rCT in receptive man
+  p1Inf_rct <- which(uCT[al[, "p1"]] == 1 & uCT.infTime[al[, "p1"]] < at &
+                       rCT[al[, "p2"]] == 0 & al[, "ins"] %in% c(1, 2))
+  p2Inf_rct <- which(uCT[al[, "p1"]] == 0 & uCT.infTime[al[, "p2"]] < at &
+                       rCT[al[, "p1"]] == 0 & al[, "ins"] %in% c(0, 2))
+  allActs_rct <- c(p1Inf_rct, p2Inf_rct)
+
+  # UAI modifier
+  uai_rct <- al[, "uai"][allActs_rct]
+  tprob_rct <- rep(rct.tprob, length(allActs_rct))
+  tprob_rct[uai_rct == 0] <- tprob_rct[uai_rct == 0] * sti.cond.rr
+
+  # Stochastic transmission
+  trans_rct <- rbinom(length(allActs_rct), 1, tprob_rct)
+
+  # Determine the newly infected partner
+  idsInf_rct <- NULL
+  if (sum(trans_rct) > 0) {
+    transAL_rct <- al[allActs_rct[trans_rct == 1],  , drop = FALSE]
+    idsInf_rct <- unique(ifelse(uCT[transAL_rct[, "p1"]] == 1,
+                                transAL_rct[, "p2"], transAL_rct[, "p1"]))
+  }
+
+  # Update attributes
+  rCT[idsInf_rct] <- 1
+  rCT.infTime[idsInf_rct] <- at
+  rCT.sympt[idsInf_rct] <- rbinom(length(idsInf_rct), 1, rct.sympt.prob)
+  rCT.timesInf[idsInf_rct] <- rCT.timesInf[idsInf_rct] + 1
+
+
+  # Urethral CT ---------------------------------------------------------
+
+  # Requires: rCT in receptive man, and no uCT in insertive man
+  p1Inf_uct <- which(rCT[al[, "p1"]] == 1 & rCT.infTime[al[, "p1"]] < at &
+                       uCT[al[, "p2"]] == 0 & al[, "ins"] %in% c(0, 2))
+  p2Inf_uct <- which(rCT[al[, "p1"]] == 0 & rCT.infTime[al[, "p2"]] < at &
+                       uCT[al[, "p1"]] == 0 & al[, "ins"] %in% c(1, 2))
+  allActs_uct <- c(p1Inf_uct, p2Inf_uct)
+
+  # UAI modifier
+  uai_uct <- al[, "uai"][allActs_uct]
+  tprob_uct <- rep(uct.tprob, length(allActs_uct))
+  tprob_uct[uai_uct == 0] <- tprob_uct[uai_uct == 0] * sti.cond.rr
+
+  # Transmission
+  trans_uct <- rbinom(length(allActs_uct), 1, tprob_uct)
+
+  # Determine the newly infected partner
+  idsInf_uct <- NULL
+  if (sum(trans_uct) > 0) {
+    transAL_uct <- al[allActs_uct[trans_uct == 1],  , drop = FALSE]
+    idsInf_uct <- unique(ifelse(uCT[transAL_uct[, "p1"]] == 1,
+                                transAL_uct[, "p2"], transAL_uct[, "p1"]))
+  }
+
+  # Update attributes
+  uCT[idsInf_uct] <- 1
+  uCT.infTime[idsInf_uct] <- at
+  uCT.sympt[idsInf_uct] <- rbinom(length(idsInf_uct), 1, uct.sympt.prob)
+  uCT.timesInf[idsInf_uct] <- uCT.timesInf[idsInf_uct] + 1
+
+
+  # Set activity cessation attribute for newly infected -----------------
+
+  # Symptomatic GC
+  GC.sympt <- which(is.na(GC.cease) & (rGC.sympt == 1 | uGC.sympt == 1))
+  idsGC.cease <- GC.sympt[which(rbinom(length(GC.sympt),
+                                       1, gc.prob.cease) == 1)]
+  GC.cease[GC.sympt] <- 0
+  GC.cease[idsGC.cease] <- 1
+
+  # Symptomatic CT
+  CT.sympt <- which(is.na(CT.cease) & (rCT.sympt == 1 | uCT.sympt == 1))
+  idsCT.cease <- CT.sympt[which(rbinom(length(CT.sympt),
+                                       1, ct.prob.cease) == 1)]
+  CT.cease[CT.sympt] <- 0
+  CT.cease[idsCT.cease] <- 1
+
+
+  # Output --------------------------------------------------------------
+
+  # attributes
+  dat$attr$rGC <- rGC
+  dat$attr$uGC <- uGC
+  dat$attr$rCT <- rCT
+  dat$attr$uCT <- uCT
+
+  dat$attr$rGC.infTime <- rGC.infTime
+  dat$attr$uGC.infTime <- uGC.infTime
+  dat$attr$rCT.infTime <- rCT.infTime
+  dat$attr$uCT.infTime <- uCT.infTime
+
+  dat$attr$rGC.timesInf <- rGC.timesInf
+  dat$attr$uGC.timesInf <- uGC.timesInf
+  dat$attr$rCT.timesInf <- rCT.timesInf
+  dat$attr$uCT.timesInf <- uCT.timesInf
+
+  dat$attr$rGC.sympt <- rGC.sympt
+  dat$attr$uGC.sympt <- uGC.sympt
+  dat$attr$rCT.sympt <- rCT.sympt
+  dat$attr$uCT.sympt <- uCT.sympt
+
+  dat$attr$GC.cease <- GC.cease
+  dat$attr$CT.cease <- CT.cease
+
+
+  # Summary stats
+  dat$epi$incid.rgc[at] <- length(idsInf_rgc)
+  dat$epi$incid.ugc[at] <- length(idsInf_ugc)
+  dat$epi$incid.gc[at] <- length(idsInf_rgc) + length(idsInf_ugc)
+  dat$epi$incid.rct[at] <- length(idsInf_rct)
+  dat$epi$incid.uct[at] <- length(idsInf_uct)
+  dat$epi$incid.ct[at] <- length(idsInf_rct) + length(idsInf_uct)
+
+  # Check all infected have all STI attributes
+  stopifnot(all(!is.na(rGC.infTime[rGC == 1])),
+            all(!is.na(rGC.sympt[rGC == 1])),
+            all(!is.na(uGC.infTime[uGC == 1])),
+            all(!is.na(uGC.sympt[uGC == 1])),
+            all(!is.na(rCT.infTime[rCT == 1])),
+            all(!is.na(rCT.sympt[rCT == 1])),
+            all(!is.na(uCT.infTime[uCT == 1])),
+            all(!is.na(uCT.sympt[uCT == 1])))
+
+  if (is.null(dat$epi$times.rgc)) {
+    dat$epi$times.rgc <- rep(NA, length(dat$epi$num))
+    dat$epi$times.ugc <- rep(NA, length(dat$epi$num))
+    dat$epi$times.rct <- rep(NA, length(dat$epi$num))
+    dat$epi$times.uct <- rep(NA, length(dat$epi$num))
+  }
+  dat$epi$times.rgc[at] <- mean(rGC.timesInf, na.rm = TRUE)
+  dat$epi$times.ugc[at] <- mean(uGC.timesInf, na.rm = TRUE)
+  dat$epi$times.rct[at] <- mean(rCT.timesInf, na.rm = TRUE)
+  dat$epi$times.uct[at] <- mean(uCT.timesInf, na.rm = TRUE)
+
+  return(dat)
+}
+
+
+#' @title STI Recovery Module
+#'
+#' @description Stochastically simulates GC/CT recovery.
+#'
+#' @inheritParams aging_msm
+#'
+#' @keywords module msm
+#'
+#' @export
+#'
+sti_recov <- function(dat, at) {
+
+  # Parameters
+  rgc.dur.asympt <- dat$param$rgc.dur.asympt
+  ugc.dur.asympt <- dat$param$ugc.dur.asympt
+  gc.dur.tx <- dat$param$gc.dur.tx
+  gc.dur.ntx <- dat$param$gc.dur.ntx
+
+  rct.dur.asympt <- dat$param$rct.dur.asympt
+  uct.dur.asympt <- dat$param$uct.dur.asympt
+  ct.dur.tx <- dat$param$ct.dur.tx
+  ct.dur.ntx <- dat$param$ct.dur.ntx
+
+
+  # GC recovery
+  idsRGC_asympt <- which(dat$attr$rGC == 1 & dat$attr$rGC.infTime < at &
+                           dat$attr$rGC.sympt == 0)
+  idsUGC_asympt <- which(dat$attr$uGC == 1 & dat$attr$uGC.infTime < at &
+                           dat$attr$uGC.sympt == 0)
+  idsRGC_tx <- which(dat$attr$rGC == 1 & dat$attr$rGC.infTime < at &
+                       dat$attr$rGC.sympt == 1 & dat$attr$rGC.tx == 1)
+  idsUGC_tx <- which(dat$attr$uGC == 1 & dat$attr$uGC.infTime < at &
+                       dat$attr$uGC.sympt == 1 & dat$attr$uGC.tx == 1)
+  idsRGC_ntx <- which(dat$attr$rGC == 1 & dat$attr$rGC.infTime < at &
+                        dat$attr$rGC.sympt == 1 & dat$attr$rGC.tx == 0)
+  idsUGC_ntx <- which(dat$attr$uGC == 1 & dat$attr$uGC.infTime < at &
+                        dat$attr$uGC.sympt == 1 & dat$attr$uGC.tx == 0)
+
+  recovRGC_asympt <- idsRGC_asympt[which(rbinom(length(idsRGC_asympt), 1,
+                                                1/rgc.dur.asympt) == 1)]
+  recovUGC_asympt <- idsUGC_asympt[which(rbinom(length(idsUGC_asympt), 1,
+                                                1/ugc.dur.asympt) == 1)]
+
+  recovRGC_tx <- idsRGC_tx[which(rbinom(length(idsRGC_tx), 1,
+                                        1/gc.dur.tx) == 1)]
+  recovUGC_tx <- idsUGC_tx[which(rbinom(length(idsUGC_tx), 1,
+                                        1/gc.dur.tx) == 1)]
+
+  if (!is.null(gc.dur.ntx)) {
+    recovRGC_ntx <- idsRGC_ntx[which(rbinom(length(idsRGC_ntx), 1,
+                                            1/gc.dur.ntx) == 1)]
+    recovUGC_ntx <- idsUGC_ntx[which(rbinom(length(idsUGC_ntx), 1,
+                                            1/gc.dur.ntx) == 1)]
+  } else {
+    recovRGC_ntx <- idsRGC_ntx[which(rbinom(length(idsRGC_ntx), 1,
+                                            1/rgc.dur.asympt) == 1)]
+    recovUGC_ntx <- idsUGC_ntx[which(rbinom(length(idsUGC_ntx), 1,
+                                            1/ugc.dur.asympt) == 1)]
+  }
+
+  recovRGC <- c(recovRGC_asympt, recovRGC_tx, recovRGC_ntx)
+  recovUGC <- c(recovUGC_asympt, recovUGC_tx, recovUGC_ntx)
+
+  dat$attr$rGC[recovRGC] <- 0
+  dat$attr$rGC.sympt[recovRGC] <- NA
+  dat$attr$rGC.infTime[recovRGC] <- NA
+  dat$attr$rGC.tx[recovRGC] <- NA
+
+  dat$attr$uGC[recovUGC] <- 0
+  dat$attr$uGC.sympt[recovUGC] <- NA
+  dat$attr$uGC.infTime[recovUGC] <- NA
+  dat$attr$uGC.tx[recovUGC] <- NA
+
+  dat$attr$GC.cease[c(recovRGC, recovUGC)] <- NA
+
+  # CT recovery
+  idsRCT_asympt <- which(dat$attr$rCT == 1 & dat$attr$rCT.infTime < at &
+                           dat$attr$rCT.sympt == 0)
+  idsUCT_asympt <- which(dat$attr$uCT == 1 & dat$attr$uCT.infTime < at &
+                           dat$attr$uCT.sympt == 0)
+  idsRCT_tx <- which(dat$attr$rCT == 1 & dat$attr$rCT.infTime < at &
+                       dat$attr$rCT.sympt == 1 & dat$attr$rCT.tx == 1)
+  idsUCT_tx <- which(dat$attr$uCT == 1 & dat$attr$uCT.infTime < at &
+                       dat$attr$uCT.sympt == 1 & dat$attr$uCT.tx == 1)
+  idsRCT_ntx <- which(dat$attr$rCT == 1 & dat$attr$rCT.infTime < at &
+                        dat$attr$rCT.sympt == 1 & dat$attr$rCT.tx == 0)
+  idsUCT_ntx <- which(dat$attr$uCT == 1 & dat$attr$uCT.infTime < at &
+                        dat$attr$uCT.sympt == 1 & dat$attr$uCT.tx == 0)
+
+  recovRCT_asympt <- idsRCT_asympt[which(rbinom(length(idsRCT_asympt),
+                                                1, 1/rct.dur.asympt) == 1)]
+  recovUCT_asympt <- idsUCT_asympt[which(rbinom(length(idsUCT_asympt),
+                                                1, 1/uct.dur.asympt) == 1)]
+
+  recovRCT_tx <- idsRCT_tx[which(rbinom(length(idsRCT_tx),
+                                        1, 1/ct.dur.tx) == 1)]
+  recovUCT_tx <- idsUCT_tx[which(rbinom(length(idsUCT_tx),
+                                        1, 1/ct.dur.tx) == 1)]
+
+  if (!is.null(ct.dur.ntx)) {
+    recovRCT_ntx <- idsRCT_ntx[which(rbinom(length(idsRCT_ntx),
+                                            1, 1/ct.dur.ntx) == 1)]
+    recovUCT_ntx <- idsUCT_ntx[which(rbinom(length(idsUCT_ntx),
+                                            1, 1/ct.dur.ntx) == 1)]
+  } else {
+    recovRCT_ntx <- idsRCT_ntx[which(rbinom(length(idsRCT_ntx),
+                                            1, 1/rct.dur.asympt) == 1)]
+    recovUCT_ntx <- idsUCT_ntx[which(rbinom(length(idsUCT_ntx),
+                                            1, 1/uct.dur.asympt) == 1)]
+  }
+
+  recovRCT <- c(recovRCT_asympt, recovRCT_tx, recovRCT_ntx)
+  recovUCT <- c(recovUCT_asympt, recovUCT_tx, recovUCT_ntx)
+
+  dat$attr$rCT[recovRCT] <- 0
+  dat$attr$rCT.sympt[recovRCT] <- NA
+  dat$attr$rCT.infTime[recovRCT] <- NA
+  dat$attr$rCT.tx[recovRCT] <- NA
+
+  dat$attr$uCT[recovUCT] <- 0
+  dat$attr$uCT.sympt[recovUCT] <- NA
+  dat$attr$uCT.infTime[recovUCT] <- NA
+  dat$attr$uCT.tx[recovUCT] <- NA
+
+  dat$attr$CT.cease[c(recovRCT, recovUCT)] <- NA
+
+  # Summary stats
+  dat$epi$recov.rgc[at] <- length(recovRGC)
+  dat$epi$recov.ugc[at] <- length(recovUGC)
+  dat$epi$recov.rct[at] <- length(recovRCT)
+  dat$epi$recov.uct[at] <- length(recovUCT)
+
+  return(dat)
+}
+
+
+#' @title STI Treatment Module
+#'
+#' @description Stochastically simulates GC/CT diagnosis and treatment.
+#'
+#' @inheritParams aging_msm
+#'
+#' @keywords module msm
+#'
+#' @export
+#'
+sti_tx <- function(dat, at) {
+
+  # Parameters
+  gc.sympt.prob.tx <- dat$param$gc.sympt.prob.tx
+  ct.sympt.prob.tx <- dat$param$ct.sympt.prob.tx
+
+  gc.asympt.prob.tx <- dat$param$gc.asympt.prob.tx
+  ct.asympt.prob.tx <- dat$param$ct.asympt.prob.tx
+
+  prep.sti.screen.int <- dat$param$prep.sti.screen.int
+  prep.sti.prob.tx <- dat$param$prep.sti.prob.tx
+
+  prep.cont.stand.tx <- dat$param$prep.continue.stand.tx
+  if (prep.cont.stand.tx == TRUE) {
+    prep.stand.tx.grp <- 0:1
+  } else {
+    prep.stand.tx.grp <- 0
+  }
+
+  # symptomatic gc treatment
+  idsRGC_tx_sympt <- which(dat$attr$rGC == 1 & dat$attr$rGC.infTime < at &
+                             dat$attr$rGC.sympt == 1 & is.na(dat$attr$rGC.tx) &
+                             dat$attr$prepStat %in% prep.stand.tx.grp)
+  idsUGC_tx_sympt <- which(dat$attr$uGC == 1 & dat$attr$uGC.infTime < at &
+                             dat$attr$uGC.sympt == 1 & is.na(dat$attr$uGC.tx) &
+                             dat$attr$prepStat %in% prep.stand.tx.grp)
+  idsGC_tx_sympt <- c(idsRGC_tx_sympt, idsUGC_tx_sympt)
+
+  txGC_sympt <- idsGC_tx_sympt[which(rbinom(length(idsGC_tx_sympt), 1,
+                                            gc.sympt.prob.tx) == 1)]
+  txRGC_sympt <- intersect(idsRGC_tx_sympt, txGC_sympt)
+  txUGC_sympt <- intersect(idsUGC_tx_sympt, txGC_sympt)
+
+  # asymptomatic gc treatment
+  idsRGC_tx_asympt <- which(dat$attr$rGC == 1 & dat$attr$rGC.infTime < at &
+                              dat$attr$rGC.sympt == 0 & is.na(dat$attr$rGC.tx) &
+                              dat$attr$prepStat %in% prep.stand.tx.grp)
+  idsUGC_tx_asympt <- which(dat$attr$uGC == 1 & dat$attr$uGC.infTime < at &
+                              dat$attr$uGC.sympt == 0 & is.na(dat$attr$uGC.tx) &
+                              dat$attr$prepStat %in% prep.stand.tx.grp)
+  idsGC_tx_asympt <- c(idsRGC_tx_asympt, idsUGC_tx_asympt)
+
+  txGC_asympt <- idsGC_tx_asympt[which(rbinom(length(idsGC_tx_asympt), 1,
+                                              gc.asympt.prob.tx) == 1)]
+  txRGC_asympt <- intersect(idsRGC_tx_asympt, txGC_asympt)
+  txUGC_asympt <- intersect(idsUGC_tx_asympt, txGC_asympt)
+
+  # all treated GC
+  txRGC <- union(txRGC_sympt, txRGC_asympt)
+  txUGC <- union(txUGC_sympt, txUGC_asympt)
+
+  idsRGC_tx <- union(idsRGC_tx_sympt, idsRGC_tx_asympt)
+  idsUGC_tx <- union(idsUGC_tx_sympt, idsUGC_tx_asympt)
+
+
+  # symptomatic ct treatment
+  idsRCT_tx_sympt <- which(dat$attr$rCT == 1 & dat$attr$rCT.infTime < at &
+                             dat$attr$rCT.sympt == 1 & is.na(dat$attr$rCT.tx) &
+                             dat$attr$prepStat %in% prep.stand.tx.grp)
+  idsUCT_tx_sympt <- which(dat$attr$uCT == 1 & dat$attr$uCT.infTime < at &
+                             dat$attr$uCT.sympt == 1 & is.na(dat$attr$uCT.tx) &
+                             dat$attr$prepStat %in% prep.stand.tx.grp)
+  idsCT_tx_sympt <- c(idsRCT_tx_sympt, idsUCT_tx_sympt)
+
+  txCT_sympt <- idsCT_tx_sympt[which(rbinom(length(idsCT_tx_sympt), 1,
+                                            ct.sympt.prob.tx) == 1)]
+  txRCT_sympt <- intersect(idsRCT_tx_sympt, txCT_sympt)
+  txUCT_sympt <- intersect(idsUCT_tx_sympt, txCT_sympt)
+
+  # asymptomatic ct treatment
+  idsRCT_tx_asympt <- which(dat$attr$rCT == 1 & dat$attr$rCT.infTime < at &
+                              dat$attr$rCT.sympt == 0 & is.na(dat$attr$rCT.tx) &
+                              dat$attr$prepStat == 0)
+  idsUCT_tx_asympt <- which(dat$attr$uCT == 1 & dat$attr$uCT.infTime < at &
+                              dat$attr$uCT.sympt == 0 & is.na(dat$attr$uCT.tx) &
+                              dat$attr$prepStat == 0)
+  idsCT_tx_asympt <- c(idsRCT_tx_asympt, idsUCT_tx_asympt)
+
+  txCT_asympt <- idsCT_tx_asympt[which(rbinom(length(idsCT_tx_asympt), 1,
+                                              ct.asympt.prob.tx) == 1)]
+  txRCT_asympt <- intersect(idsRCT_tx_asympt, txCT_asympt)
+  txUCT_asympt <- intersect(idsUCT_tx_asympt, txCT_asympt)
+
+  # all treated CT
+  txRCT <- union(txRCT_sympt, txRCT_asympt)
+  txUCT <- union(txUCT_sympt, txUCT_asympt)
+
+  idsRCT_tx <- union(idsRCT_tx_sympt, idsRCT_tx_asympt)
+  idsUCT_tx <- union(idsUCT_tx_sympt, idsUCT_tx_asympt)
+
+  # Interval-based treatment for MSM on PrEP
+  idsSTI_screen <- which(dat$attr$prepStartTime == at |
+                           (at - dat$attr$prepLastStiScreen >= prep.sti.screen.int))
+
+  dat$attr$prepLastStiScreen[idsSTI_screen] <- at
+
+  idsRGC_prep_tx <- intersect(idsSTI_screen,
+                              which(dat$attr$rGC == 1 &
+                                      dat$attr$rGC.infTime < at &
+                                      is.na(dat$attr$rGC.tx)))
+  idsUGC_prep_tx <- intersect(idsSTI_screen,
+                              which(dat$attr$uGC == 1 &
+                                      dat$attr$uGC.infTime < at &
+                                      is.na(dat$attr$uGC.tx)))
+  idsRCT_prep_tx <- intersect(idsSTI_screen,
+                              which(dat$attr$rCT == 1 &
+                                      dat$attr$rCT.infTime < at &
+                                      is.na(dat$attr$rCT.tx)))
+  idsUCT_prep_tx <- intersect(idsSTI_screen,
+                              which(dat$attr$uCT == 1 & dat$attr$uCT.infTime < at &
+                                      is.na(dat$attr$uCT.tx)))
+
+  txRGC_prep <- idsRGC_prep_tx[which(rbinom(length(idsRGC_prep_tx), 1,
+                                            prep.sti.prob.tx) == 1)]
+  txUGC_prep <- idsUGC_prep_tx[which(rbinom(length(idsUGC_prep_tx), 1,
+                                            prep.sti.prob.tx) == 1)]
+  txRCT_prep <- idsRCT_prep_tx[which(rbinom(length(idsRCT_prep_tx), 1,
+                                            prep.sti.prob.tx) == 1)]
+  txUCT_prep <- idsUCT_prep_tx[which(rbinom(length(idsUCT_prep_tx), 1,
+                                            prep.sti.prob.tx) == 1)]
+
+  txRGC <- c(txRGC, txRGC_prep)
+  txUGC <- c(txUGC, txUGC_prep)
+  txRCT <- c(txRCT, txRCT_prep)
+  txUCT <- c(txUCT, txUCT_prep)
+
+
+  # update attributes
+  dat$attr$rGC.tx[c(idsRGC_tx, idsRGC_prep_tx)] <- 0
+  dat$attr$rGC.tx[txRGC] <- 1
+
+  dat$attr$uGC.tx[c(idsUGC_tx, idsUGC_prep_tx)] <- 0
+  dat$attr$uGC.tx[txUGC] <- 1
+
+  dat$attr$rCT.tx[c(idsRCT_tx, idsRCT_prep_tx)] <- 0
+  dat$attr$rCT.tx[txRCT] <- 1
+
+  dat$attr$uCT.tx[c(idsUCT_tx, idsUCT_prep_tx)] <- 0
+  dat$attr$uCT.tx[txUCT] <- 1
+
+  # add tx at other site
+  dat$attr$rGC.tx[which(dat$attr$uGC.tx == 1 & dat$attr$rGC == 1)] <- 1
+  dat$attr$uGC.tx[which(dat$attr$rGC.tx == 1 & dat$attr$rGC == 1)] <- 1
+
+  dat$attr$rCT.tx[which(dat$attr$uCT.tx == 1 & dat$attr$rCT == 1)] <- 1
+  dat$attr$uCT.tx[which(dat$attr$rCT.tx == 1 & dat$attr$rCT == 1)] <- 1
+
+  # summary stats
+  if (is.null(dat$epi$txGC)) {
+    dat$epi$txGC <- rep(NA, length(dat$epi$num))
+    dat$epi$txCT <- rep(NA, length(dat$epi$num))
+  }
+  dat$epi$txGC[at] <- length(txRGC) + length(txUGC)
+  dat$epi$txCT[at] <- length(txRCT) + length(txUCT)
+
+  return(dat)
+}
diff --git a/R/mod.trans.R b/R/mod.trans.R
index 99e2ffaf..bc99eae5 100644
--- a/R/mod.trans.R
+++ b/R/mod.trans.R
@@ -32,7 +32,7 @@
 #'
 #' @export
 #'
-trans_msm <- function(dat, at){
+trans_msm <- function(dat, at) {
 
   # Variables -----------------------------------------------------------
 
@@ -41,10 +41,13 @@ trans_msm <- function(dat, at){
   stage <- dat$attr$stage
   ccr5 <- dat$attr$ccr5
   circ <- dat$attr$circ
-  diag.status <- dat$attr$diag.status
-  tx.status <- dat$attr$tx.status
+  status <- dat$attr$status
   prepStat <- dat$attr$prepStat
   prepClass <- dat$attr$prepClass
+  rGC <- dat$attr$rGC
+  uGC <- dat$attr$uGC
+  rCT <- dat$attr$rCT
+  uCT <- dat$attr$uCT
 
   # Parameters
   URAI.prob <- dat$param$URAI.prob
@@ -54,9 +57,16 @@ trans_msm <- function(dat, at){
   circ.rr <- dat$param$circ.rr
   ccr5.heteroz.rr <- dat$param$ccr5.heteroz.rr
   prep.hr <- dat$param$prep.class.hr
+  hiv.ugc.rr <- dat$param$hiv.ugc.rr
+  hiv.uct.rr <- dat$param$hiv.uct.rr
+  hiv.rgc.rr <- dat$param$hiv.rgc.rr
+  hiv.rct.rr <- dat$param$hiv.rct.rr
+  hiv.dual.rr <- dat$param$hiv.dual.rr
+
 
   # Data
-  dal <- dat$temp$dal
+  al <- dat$temp$al
+  dal <- al[which(status[al[, 1]] == 1 & status[al[, 2]] == 0), ]
   dal <- dal[sample(1:nrow(dal)), ]
   ncols <- dim(dal)[2]
 
@@ -64,12 +74,10 @@ trans_msm <- function(dat, at){
     return(dat)
   }
 
-  ## Reorder by role: ins on the left, rec on the right,
-  ##                  with flippers represented twice
+  ## Reorder by role: ins on the left, rec on the right, flippers represented twice
   disc.ip <- dal[dal[, "ins"] %in% 1:2, ]
   disc.rp <- dal[dal[, "ins"] %in% c(0, 2), c(2:1, 3:ncols)]
-  colnames(disc.ip)[1:2] <- c("i", "r")
-  colnames(disc.rp)[1:2] <- c("i", "r")
+  colnames(disc.ip)[1:2] <- colnames(disc.rp)[1:2] <- c("ins", "rec")
 
 
   # PATP: Insertive Man Infected (Col 1) --------------------------------
@@ -82,6 +90,8 @@ trans_msm <- function(dat, at){
   ip.ccr5 <- ccr5[disc.ip[, 2]]
   ip.prep <- prepStat[disc.ip[, 2]]
   ip.prepcl <- prepClass[disc.ip[, 2]]
+  ip.rGC <- rGC[disc.ip[, 2]]
+  ip.rCT <- rCT[disc.ip[, 2]]
 
   # Base TP from VL
   ip.tprob <- URAI.prob * 2.45^(ip.vl - 4.5)
@@ -104,10 +114,26 @@ trans_msm <- function(dat, at){
   }
 
   # Acute-stage multipliers
-  isAcute <- which(ip.stage %in% c(1, 2))
+  isAcute <- which(ip.stage %in% 1:2)
   ip.tlo[isAcute] <- ip.tlo[isAcute] + log(acute.rr)
 
-  # Retransformation to probability
+  ## Multiplier for STI
+  is.rGC <- which(ip.rGC == 1)
+
+  is.rCT <- which(ip.rCT == 1)
+
+  is.rect.dual <- intersect(is.rGC, is.rCT)
+
+  is.rGC.sing <- setdiff(is.rGC, is.rect.dual)
+  is.rCT.sing <- setdiff(is.rCT, is.rect.dual)
+
+  ip.tlo[is.rGC.sing] <- ip.tlo[is.rGC.sing] + log(hiv.rgc.rr)
+  ip.tlo[is.rCT.sing] <- ip.tlo[is.rCT.sing] + log(hiv.rct.rr)
+
+  ip.tlo[is.rect.dual] <- ip.tlo[is.rect.dual] +
+    max(log(hiv.rgc.rr), log(hiv.rct.rr)) +
+    min(log(hiv.rgc.rr), log(hiv.rct.rr)) * hiv.dual.rr
+
   ip.tprob <- plogis(ip.tlo)
   stopifnot(ip.tprob >= 0, ip.tprob <= 1)
 
@@ -123,6 +149,8 @@ trans_msm <- function(dat, at){
   rp.ccr5 <- ccr5[disc.rp[, 1]]
   rp.prep <- prepStat[disc.rp[, 1]]
   rp.prepcl <- prepClass[disc.rp[, 1]]
+  rp.uGC <- uGC[disc.rp[, 1]]
+  rp.uCT <- uCT[disc.rp[, 1]]
 
   # Base TP from VL
   rp.tprob <- UIAI.prob * 2.45^(rp.vl - 4.5)
@@ -148,13 +176,31 @@ trans_msm <- function(dat, at){
   }
 
   # Acute-stage multipliers
-  isAcute <- which(rp.stage %in% c(1, 2))
+  isAcute <- which(rp.stage %in% 1:2)
   rp.tlo[isAcute] <- rp.tlo[isAcute] + log(acute.rr)
 
+  ## Multiplier for STI
+  is.uGC <- which(rp.uGC == 1)
+
+  is.uCT <- which(rp.uCT == 1)
+
+  is.ureth.dual <- intersect(is.uGC, is.uCT)
+
+  is.uGC.sing <- setdiff(is.uGC, is.ureth.dual)
+  is.uCT.sing <- setdiff(is.uCT, is.ureth.dual)
+
+  rp.tlo[is.uGC.sing] <- rp.tlo[is.uGC.sing] + log(hiv.ugc.rr)
+  rp.tlo[is.uCT.sing] <- rp.tlo[is.uCT.sing] + log(hiv.uct.rr)
+
+  rp.tlo[is.ureth.dual] <- rp.tlo[is.ureth.dual] +
+    max(log(hiv.ugc.rr), log(hiv.uct.rr)) +
+    min(log(hiv.ugc.rr), log(hiv.uct.rr)) * hiv.dual.rr
+
   # Retransformation to probability
   rp.tprob <- plogis(rp.tlo)
   stopifnot(rp.tprob >= 0, rp.tprob <= 1)
 
+
   # Transmission --------------------------------------------------------
 
   ## Bernoulli transmission events
@@ -164,23 +210,18 @@ trans_msm <- function(dat, at){
 
   # Output --------------------------------------------------------------
 
+
   # Update attributes
 
-  infected <- infector <- inf.type <- NULL
+  infected <- inf.type <- NULL
   if (sum(trans.ip, trans.rp) > 0) {
 
     infected <- c(disc.ip[trans.ip == 1, 2],
                   disc.rp[trans.rp == 1, 1])
-    infector <- c(disc.ip[trans.ip == 1, 1],
-                  disc.rp[trans.rp == 1, 2])
     inf.role <- c(rep(0, sum(trans.ip)), rep(1, sum(trans.rp)))
     inf.type <- c(disc.ip[trans.ip == 1, "ptype"],
                   disc.rp[trans.rp == 1, "ptype"])
 
-    inf.stage <- stage[infector]
-    inf.diag <- diag.status[infector]
-    inf.tx <- tx.status[infector]
-
     dat$attr$status[infected] <- 1
     dat$attr$inf.time[infected] <- at
     dat$attr$vl[infected] <- 0
@@ -189,12 +230,8 @@ trans_msm <- function(dat, at){
     dat$attr$diag.status[infected] <- 0
     dat$attr$tx.status[infected] <- 0
 
-    dat$attr$infector[infected] <- infector
     dat$attr$inf.role[infected] <- inf.role
     dat$attr$inf.type[infected] <- inf.type
-    dat$attr$inf.diag[infected] <- inf.diag
-    dat$attr$inf.tx[infected] <- inf.tx
-    dat$attr$inf.stage[infected] <- inf.stage
 
     dat$attr$cum.time.on.tx[infected] <- 0
     dat$attr$cum.time.off.tx[infected] <- 0
@@ -203,6 +240,9 @@ trans_msm <- function(dat, at){
   # Summary Output
   dat$epi$incid[at] <- length(infected)
 
+  dat$epi$trans.main[at] <- sum(inf.type == 1)
+  dat$epi$trans.casl[at] <- sum(inf.type == 2)
+  dat$epi$trans.inst[at] <- sum(inf.type == 3)
 
   return(dat)
 }
diff --git a/R/mod.verbose.R b/R/mod.verbose.R
index 99f502e2..d76925b0 100644
--- a/R/mod.verbose.R
+++ b/R/mod.verbose.R
@@ -19,7 +19,7 @@
 #' exist.
 #'
 #' @keywords module msm
-#' 
+#'
 #' @export
 #'
 verbose_msm <- function(x, type, s, at) {
@@ -44,8 +44,6 @@ verbose_msm <- function(x, type, s, at) {
                 "\nStep: ", at, " (", round(at/x$control$nsteps, 2), ")",
                 "\nPop Size: ", x$epi$num[at],
                 "\nTot Prev: ", round(x$epi$i.num[at] / x$epi$num[at], 3),
-                "\nWht Prev: ", round(x$epi$i.num.W[at] / x$epi$num.W[at], 3),
-                "\nBlk Prev: ", round(x$epi$i.num.B[at] / x$epi$num.B[at], 3),
                 "\n\n", sep = "", file = fn)
           }
         }
@@ -56,8 +54,10 @@ verbose_msm <- function(x, type, s, at) {
           nsteps <- x$control$nsteps
           time.unit <- x$param$time.unit
           prev <- round(x$epi$i.prev[at], 3)
-          prev.B <- round(x$epi$i.prev.B[at], 3)
-          prev.W <- round(x$epi$i.prev.W[at], 3)
+          prev.rgc <- round(x$epi$prev.rgc[at], 3)
+          prev.ugc <- round(x$epi$prev.ugc[at], 3)
+          prev.rct <- round(x$epi$prev.rct[at], 3)
+          prev.uct <- round(x$epi$prev.uct[at], 3)
 
           cat("\014")
           cat("\nEpidemic Simulation")
@@ -68,14 +68,15 @@ verbose_msm <- function(x, type, s, at) {
               round((nsteps * time.unit) / 365, 1), sep = "")
           cat("\n------------------------------")
           cat("\nTotal Pop Size:", x$epi$num[at])
-          cat("\nBlack Pop Size:", x$epi$num.B[at])
-          cat("\nWhite Pop Size:", x$epi$num.W[at])
           cat("\n------------------------------")
-          cat("\nCurr Incidence:", x$epi$incid[at])
-          cat("\nCuml Incidence:", sum(x$epi$incid, na.rm = TRUE))
-          cat("\nTotal Prevalence: ", x$epi$i.num[at], " (", prev, ")", sep = "")
-          cat("\nBlack Prevalence: ", x$epi$i.num.B[at], " (", prev.B, ")", sep = "")
-          cat("\nWhite Prevalence: ", x$epi$i.num.W[at], " (", prev.W, ")", sep = "")
+          cat("\nHIV Curr Incidence:", x$epi$incid[at])
+          cat("\nHIV Cuml Incidence:", sum(x$epi$incid, na.rm = TRUE))
+          cat("\nHIV Prevalence: ", x$epi$i.num[at], " (", prev, ")", sep = "")
+          cat("\n------------------------------")
+          cat("\nrGC Prevalence: ", prev.rgc, sep = "")
+          cat("\nuGC Prevalence: ", prev.ugc, sep = "")
+          cat("\nrCT Prevalence: ", prev.rct, sep = "")
+          cat("\nuCT Prevalence: ", prev.uct, sep = "")
           cat("\n==============================")
 
         }
@@ -107,7 +108,7 @@ verbose_msm <- function(x, type, s, at) {
 #' exist.
 #'
 #' @keywords module het
-#' 
+#'
 #' @export
 #'
 verbose_het <- function(x, type, s, at) {
diff --git a/R/params.R b/R/params.R
index b656aa10..beb244cf 100644
--- a/R/params.R
+++ b/R/params.R
@@ -205,6 +205,7 @@
 #' @param vv.iev.WW.prob Probability that in a white-white partnership of
 #'        two versatile men, they will engage in intra-event versatility
 #'        ("flipping") given that they're having AI.
+#'
 #' @param prep.start Time step at which the PrEP intervention should start.
 #' @param prep.elig.model Modeling approach for determining who is eligible for
 #'        PrEP. Current options are limited to: \code{"all"} for all persons who
@@ -227,6 +228,63 @@
 #'        in days.
 #' @param prep.risk.reassess If \code{TRUE}, reassess eligibility for PrEP at
 #'        each testing visit.
+#'
+#' @param rcomp.prob Level of risk compensation from 0 to 1, where 0 is no risk
+#'        compensation, 0.5 is a 50% reduction in the probability of condom use
+#'        per act, and 1 is a complete cessation of condom use following PrEP
+#'        initiation.
+#' @param rcomp.adh.groups PrEP adherence groups for whom risk compensation
+#'        occurs, as a vector with values 0, 1, 2, 3 corresponding to non-adherent,
+#'        low adherence, medium adherence, and high adherence to PrEP.
+#' @param rcomp.main.only Logical, if risk compensation is limited to main
+#'        partnerships only, versus all partnerships.
+#' @param rcomp.discl.only Logical, if risk compensation is limited known-discordant
+#'        partnerships only, versus all partnerships.
+#'
+#' @param rgc.tprob Probability of rectal gonorrhea infection per act.
+#' @param ugc.tprob Probability of urethral gonorrhea infection per act.
+#' @param rct.tprob Probability of rectal chylamdia infection per act.
+#' @param uct.tprob Probability of urethral chylamdia infection per act.
+#' @param rgc.sympt.prob Probability of symptoms given infection with rectal
+#'        gonorrhea.
+#' @param ugc.sympt.prob Probability of symptoms given infection with urethral
+#'        gonorrhea.
+#' @param rct.sympt.prob Probability of symptoms given infection with rectal
+#'        chylamdia.
+#' @param uct.sympt.prob Probability of symptoms given infection with urethral
+#'        chylamdia.
+#' @param rgc.dur.asympt Average duration in weeks of asymptomatic rectal gonorrhea.
+#' @param ugc.dur.asympt Average duration in weeks of asymptomatic urethral gonorrhea.
+#' @param gc.dur.tx Average duration in weeks of treated gonorrhea (both sites).
+#' @param gc.dur.ntx Average duration in weeks of untreated, symptomatic gonorrhea (both sites).
+#'        If \code{NULL}, uses site-specific durations for asymptomatic infections.
+#' @param rct.dur.asympt Average in weeks duration of asymptomatic rectal chylamdia.
+#' @param uct.dur.asympt Average in weeks duration of asymptomatic urethral chylamdia.
+#' @param ct.dur.tx Average in weeks duration of treated chylamdia (both sites).
+#' @param ct.dur.ntx Average in weeks duration of untreated, symptomatic chylamdia (both sites).
+#'        If \code{NULL}, uses site-specific durations for asymptomatic infections.
+#' @param gc.prob.cease Probability of ceasing sexual activity during symptomatic
+#'        infection with gonorrhea.
+#' @param ct.prob.cease Probability of ceasing sexual activity during symptomatic
+#'        infection with chylamdia.
+#' @param gc.sympt.prob.tx Probability of treatment for symptomatic gonorrhea.
+#' @param ct.sympt.prob.tx Probability of treatment for symptomatic chylamdia.
+#' @param gc.asympt.prob.tx Probability of treatment for asymptomatic gonorrhea.
+#' @param ct.asympt.prob.tx Probability of treatment for asymptomatic chylamdia.
+#' @param prep.sti.screen.int Interval in days between STI screening at PrEP visits.
+#' @param prep.sti.prob.tx Probability of treatment given positive screening during
+#'        PrEP visit.
+#' @param prep.continue.stand.tx Logical, if \code{TRUE} will continue standard
+#'        STI treatment of symptomatic cases even after PrEP initiation.
+#' @param sti.cond.rr Relative risk of STI infection (in either direction) given
+#'        a condom used by the insertive partner.
+#' @param hiv.rgc.rr Relative risk of HIV infection given current rectal gonorrhea.
+#' @param hiv.ugc.rr Relative risk of HIV infection given current urethral gonorrhea.
+#' @param hiv.rct.rr Relative risk of HIV infection given current rectal chylamdia.
+#' @param hiv.uct.rr Relative risk of HIV infection given current urethral chylamdia.
+#' @param hiv.dual.rr Additive proportional risk, from 0 to 1, for HIV infection
+#'        given dual infection with both gonorrhea and chylamdia.
+#'
 #' @param ... Additional arguments passed to the function.
 #'
 #' @return
@@ -236,6 +294,7 @@
 #' @keywords msm
 #'
 #' @export
+#'
 param_msm <- function(nwstats,
                       race.method = 1,
                       last.neg.test.B.int = 301,
@@ -350,6 +409,51 @@ param_msm <- function(nwstats,
                       prep.tst.int = 90,
                       prep.risk.int = 182,
                       prep.risk.reassess = TRUE,
+
+                      rcomp.prob = 0,
+                      rcomp.adh.groups = 0:3,
+                      rcomp.main.only = FALSE,
+                      rcomp.discl.only = FALSE,
+
+                      rgc.tprob = 0.45717564,
+                      ugc.tprob = 0.18828959,
+                      rct.tprob = 0.40582277,
+                      uct.tprob = 0.16310650,
+
+                      rgc.sympt.prob = 0.07427077,
+                      ugc.sympt.prob = 0.86191181,
+                      rct.sympt.prob = 0.08108045,
+                      uct.sympt.prob = 0.91586663,
+
+                      rgc.dur.asympt = 42.82357437,
+                      ugc.dur.asympt = 33.43336473,
+                      gc.dur.tx = 2,
+                      gc.dur.ntx = NULL,
+
+                      rct.dur.asympt = 54.80423283,
+                      uct.dur.asympt = 46.20741625,
+                      ct.dur.tx = 2,
+                      ct.dur.ntx = NULL,
+
+                      gc.prob.cease = 0.03344859,
+                      ct.prob.cease = 0.03344859,
+
+                      gc.sympt.prob.tx = 0.90,
+                      ct.sympt.prob.tx = 0.85,
+                      gc.asympt.prob.tx = 0,
+                      ct.asympt.prob.tx = 0,
+
+                      prep.sti.screen.int = 182,
+                      prep.sti.prob.tx = 1,
+                      prep.continue.stand.tx = TRUE,
+
+                      sti.cond.rr = 0.3,
+
+                      hiv.rgc.rr = 2.83399215,
+                      hiv.ugc.rr = 1.56832041,
+                      hiv.rct.rr = 2.83399215,
+                      hiv.uct.rr = 1.56832041,
+                      hiv.dual.rr = 0,
                       ...) {
 
   p <- get_args(formal.args = formals(sys.function()),
@@ -519,8 +623,13 @@ init_msm <- function(nwstats,
 #' @param riskhist.FUN Module function to calculate risk history for uninfected
 #'        persons in the population.
 #' @param position.FUN Module function to simulate sexual position within acts.
-#' @param trans.FUN Module function to stochastically simulate disease transmission
+#' @param trans.FUN Module function to stochastically simulate HIV transmission
 #'        over acts given individual and dyadic attributes.
+#' @param stitrans.FUN Module function to simulate GC/CT transmission over current
+#'        edgelist.
+#' @param stirecov.FUN Module function to simulate recovery from GC/CT, heterogeneous
+#'        by disease, site, symptoms, and treatment status.
+#' @param stitx.FUN Module function to simulate treatment of GC/CT.
 #' @param prev.FUN Module function to calculate prevalence summary statistics.
 #' @param verbose.FUN Module function to print model progress to the console or
 #'        external text files.
@@ -564,6 +673,9 @@ control_msm <- function(simno = 1,
                         riskhist.FUN = riskhist_msm,
                         position.FUN = position_msm,
                         trans.FUN = trans_msm,
+                        stitrans.FUN = sti_trans,
+                        stirecov.FUN = sti_recov,
+                        stitx.FUN = sti_tx,
                         prev.FUN = prevalence_msm,
                         verbose.FUN = verbose_msm,
                         save.nwstats = FALSE,
diff --git a/man/control_msm.Rd b/man/control_msm.Rd
index 88cc810a..3a033034 100644
--- a/man/control_msm.Rd
+++ b/man/control_msm.Rd
@@ -12,7 +12,8 @@ control_msm(simno = 1, nsims = 1, ncores = 1, nsteps = 100, start = 1,
   resim_nets.FUN = simnet_msm, disclose.FUN = disclose_msm,
   acts.FUN = acts_msm, condoms.FUN = condoms_msm,
   riskhist.FUN = riskhist_msm, position.FUN = position_msm,
-  trans.FUN = trans_msm, prev.FUN = prevalence_msm,
+  trans.FUN = trans_msm, stitrans.FUN = sti_trans,
+  stirecov.FUN = sti_recov, stitx.FUN = sti_tx, prev.FUN = prevalence_msm,
   verbose.FUN = verbose_msm, save.nwstats = FALSE, save.other = c("attr",
   "temp", "el", "p"), verbose = TRUE, verbose.int = 1, ...)
 }
@@ -70,9 +71,17 @@ persons in the population.}
 
 \item{position.FUN}{Module function to simulate sexual position within acts.}
 
-\item{trans.FUN}{Module function to stochastically simulate disease transmission
+\item{trans.FUN}{Module function to stochastically simulate HIV transmission
 over acts given individual and dyadic attributes.}
 
+\item{stitrans.FUN}{Module function to simulate GC/CT transmission over current
+edgelist.}
+
+\item{stirecov.FUN}{Module function to simulate recovery from GC/CT, heterogeneous
+by disease, site, symptoms, and treatment status.}
+
+\item{stitx.FUN}{Module function to simulate treatment of GC/CT.}
+
 \item{prev.FUN}{Module function to calculate prevalence summary statistics.}
 
 \item{verbose.FUN}{Module function to print model progress to the console or
diff --git a/man/param_msm.Rd b/man/param_msm.Rd
index 611c7a3f..39ac31eb 100644
--- a/man/param_msm.Rd
+++ b/man/param_msm.Rd
@@ -47,7 +47,21 @@ param_msm(nwstats, race.method = 1, last.neg.test.B.int = 301,
   prep.elig.model = "base", prep.class.prob = c(0.211, 0.07, 0.1, 0.619),
   prep.class.hr = c(1, 0.69, 0.19, 0.05), prep.coverage = 0,
   prep.cov.method = "curr", prep.cov.rate = 1, prep.tst.int = 90,
-  prep.risk.int = 182, prep.risk.reassess = TRUE, ...)
+  prep.risk.int = 182, prep.risk.reassess = TRUE, rcomp.prob = 0,
+  rcomp.adh.groups = 0:3, rcomp.main.only = FALSE,
+  rcomp.discl.only = FALSE, rgc.tprob = 0.45717564,
+  ugc.tprob = 0.18828959, rct.tprob = 0.40582277, uct.tprob = 0.1631065,
+  rgc.sympt.prob = 0.07427077, ugc.sympt.prob = 0.86191181,
+  rct.sympt.prob = 0.08108045, uct.sympt.prob = 0.91586663,
+  rgc.dur.asympt = 42.82357437, ugc.dur.asympt = 33.43336473,
+  gc.dur.tx = 2, gc.dur.ntx = NULL, rct.dur.asympt = 54.80423283,
+  uct.dur.asympt = 46.20741625, ct.dur.tx = 2, ct.dur.ntx = NULL,
+  gc.prob.cease = 0.03344859, ct.prob.cease = 0.03344859,
+  gc.sympt.prob.tx = 0.9, ct.sympt.prob.tx = 0.85, gc.asympt.prob.tx = 0,
+  ct.asympt.prob.tx = 0, prep.sti.screen.int = 182, prep.sti.prob.tx = 1,
+  prep.continue.stand.tx = TRUE, sti.cond.rr = 0.3,
+  hiv.rgc.rr = 2.83399215, hiv.ugc.rr = 1.56832041,
+  hiv.rct.rr = 2.83399215, hiv.uct.rr = 1.56832041, hiv.dual.rr = 0, ...)
 }
 \arguments{
 \item{nwstats}{Target statistics for the network model. An object of class
@@ -373,6 +387,95 @@ in days.}
 \item{prep.risk.reassess}{If \code{TRUE}, reassess eligibility for PrEP at
 each testing visit.}
 
+\item{rcomp.prob}{Level of risk compensation from 0 to 1, where 0 is no risk
+compensation, 0.5 is a 50% reduction in the probability of condom use
+per act, and 1 is a complete cessation of condom use following PrEP
+initiation.}
+
+\item{rcomp.adh.groups}{PrEP adherence groups for whom risk compensation
+occurs, as a vector with values 0, 1, 2, 3 corresponding to non-adherent,
+low adherence, medium adherence, and high adherence to PrEP.}
+
+\item{rcomp.main.only}{Logical, if risk compensation is limited to main
+partnerships only, versus all partnerships.}
+
+\item{rcomp.discl.only}{Logical, if risk compensation is limited known-discordant
+partnerships only, versus all partnerships.}
+
+\item{rgc.tprob}{Probability of rectal gonorrhea infection per act.}
+
+\item{ugc.tprob}{Probability of urethral gonorrhea infection per act.}
+
+\item{rct.tprob}{Probability of rectal chylamdia infection per act.}
+
+\item{uct.tprob}{Probability of urethral chylamdia infection per act.}
+
+\item{rgc.sympt.prob}{Probability of symptoms given infection with rectal
+gonorrhea.}
+
+\item{ugc.sympt.prob}{Probability of symptoms given infection with urethral
+gonorrhea.}
+
+\item{rct.sympt.prob}{Probability of symptoms given infection with rectal
+chylamdia.}
+
+\item{uct.sympt.prob}{Probability of symptoms given infection with urethral
+chylamdia.}
+
+\item{rgc.dur.asympt}{Average duration in weeks of asymptomatic rectal gonorrhea.}
+
+\item{ugc.dur.asympt}{Average duration in weeks of asymptomatic urethral gonorrhea.}
+
+\item{gc.dur.tx}{Average duration in weeks of treated gonorrhea (both sites).}
+
+\item{gc.dur.ntx}{Average duration in weeks of untreated, symptomatic gonorrhea (both sites).
+If \code{NULL}, uses site-specific durations for asymptomatic infections.}
+
+\item{rct.dur.asympt}{Average in weeks duration of asymptomatic rectal chylamdia.}
+
+\item{uct.dur.asympt}{Average in weeks duration of asymptomatic urethral chylamdia.}
+
+\item{ct.dur.tx}{Average in weeks duration of treated chylamdia (both sites).}
+
+\item{ct.dur.ntx}{Average in weeks duration of untreated, symptomatic chylamdia (both sites).
+If \code{NULL}, uses site-specific durations for asymptomatic infections.}
+
+\item{gc.prob.cease}{Probability of ceasing sexual activity during symptomatic
+infection with gonorrhea.}
+
+\item{ct.prob.cease}{Probability of ceasing sexual activity during symptomatic
+infection with chylamdia.}
+
+\item{gc.sympt.prob.tx}{Probability of treatment for symptomatic gonorrhea.}
+
+\item{ct.sympt.prob.tx}{Probability of treatment for symptomatic chylamdia.}
+
+\item{gc.asympt.prob.tx}{Probability of treatment for asymptomatic gonorrhea.}
+
+\item{ct.asympt.prob.tx}{Probability of treatment for asymptomatic chylamdia.}
+
+\item{prep.sti.screen.int}{Interval in days between STI screening at PrEP visits.}
+
+\item{prep.sti.prob.tx}{Probability of treatment given positive screening during
+PrEP visit.}
+
+\item{prep.continue.stand.tx}{Logical, if \code{TRUE} will continue standard
+STI treatment of symptomatic cases even after PrEP initiation.}
+
+\item{sti.cond.rr}{Relative risk of STI infection (in either direction) given
+a condom used by the insertive partner.}
+
+\item{hiv.rgc.rr}{Relative risk of HIV infection given current rectal gonorrhea.}
+
+\item{hiv.ugc.rr}{Relative risk of HIV infection given current urethral gonorrhea.}
+
+\item{hiv.rct.rr}{Relative risk of HIV infection given current rectal chylamdia.}
+
+\item{hiv.uct.rr}{Relative risk of HIV infection given current urethral chylamdia.}
+
+\item{hiv.dual.rr}{Additive proportional risk, from 0 to 1, for HIV infection
+given dual infection with both gonorrhea and chylamdia.}
+
 \item{...}{Additional arguments passed to the function.}
 }
 \value{
diff --git a/man/sti_recov.Rd b/man/sti_recov.Rd
new file mode 100644
index 00000000..03707965
--- /dev/null
+++ b/man/sti_recov.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mod.sti.R
+\name{sti_recov}
+\alias{sti_recov}
+\title{STI Recovery Module}
+\usage{
+sti_recov(dat, at)
+}
+\arguments{
+\item{dat}{Master data list object of class \code{dat} containing networks,
+individual-level attributes, and summary statistics.}
+
+\item{at}{Current time step.}
+}
+\description{
+Stochastically simulates GC/CT recovery.
+}
+\keyword{module}
+\keyword{msm}
+
diff --git a/man/sti_trans.Rd b/man/sti_trans.Rd
new file mode 100644
index 00000000..b2ff0bff
--- /dev/null
+++ b/man/sti_trans.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mod.sti.R
+\name{sti_trans}
+\alias{sti_trans}
+\title{STI Transmission Module}
+\usage{
+sti_trans(dat, at)
+}
+\arguments{
+\item{dat}{Master data list object of class \code{dat} containing networks,
+individual-level attributes, and summary statistics.}
+
+\item{at}{Current time step.}
+}
+\description{
+Stochastically simulates GC/CT transmission given the current
+             state of the edgelist.
+}
+\keyword{module}
+\keyword{msm}
+
diff --git a/man/sti_tx.Rd b/man/sti_tx.Rd
new file mode 100644
index 00000000..f5b88ebd
--- /dev/null
+++ b/man/sti_tx.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mod.sti.R
+\name{sti_tx}
+\alias{sti_tx}
+\title{STI Treatment Module}
+\usage{
+sti_tx(dat, at)
+}
+\arguments{
+\item{dat}{Master data list object of class \code{dat} containing networks,
+individual-level attributes, and summary statistics.}
+
+\item{at}{Current time step.}
+}
+\description{
+Stochastically simulates GC/CT diagnosis and treatment.
+}
+\keyword{module}
+\keyword{msm}
+

From ac539b429b07249e93f05320465c9b5859a5c657 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Tue, 9 Aug 2016 10:05:26 -0400
Subject: [PATCH 02/40] Fix typos

---
 R/params.R       | 40 ++++++++++++++++++++++++----------------
 man/init_msm.Rd  | 11 ++++++++++-
 man/param_msm.Rd | 28 ++++++++++++++--------------
 3 files changed, 48 insertions(+), 31 deletions(-)

diff --git a/R/params.R b/R/params.R
index beb244cf..0984e1d7 100644
--- a/R/params.R
+++ b/R/params.R
@@ -243,34 +243,34 @@
 #'
 #' @param rgc.tprob Probability of rectal gonorrhea infection per act.
 #' @param ugc.tprob Probability of urethral gonorrhea infection per act.
-#' @param rct.tprob Probability of rectal chylamdia infection per act.
-#' @param uct.tprob Probability of urethral chylamdia infection per act.
+#' @param rct.tprob Probability of rectal chlamydia infection per act.
+#' @param uct.tprob Probability of urethral chlamydia infection per act.
 #' @param rgc.sympt.prob Probability of symptoms given infection with rectal
 #'        gonorrhea.
 #' @param ugc.sympt.prob Probability of symptoms given infection with urethral
 #'        gonorrhea.
 #' @param rct.sympt.prob Probability of symptoms given infection with rectal
-#'        chylamdia.
+#'        chlamydia.
 #' @param uct.sympt.prob Probability of symptoms given infection with urethral
-#'        chylamdia.
+#'        chlamydia.
 #' @param rgc.dur.asympt Average duration in weeks of asymptomatic rectal gonorrhea.
 #' @param ugc.dur.asympt Average duration in weeks of asymptomatic urethral gonorrhea.
 #' @param gc.dur.tx Average duration in weeks of treated gonorrhea (both sites).
 #' @param gc.dur.ntx Average duration in weeks of untreated, symptomatic gonorrhea (both sites).
 #'        If \code{NULL}, uses site-specific durations for asymptomatic infections.
-#' @param rct.dur.asympt Average in weeks duration of asymptomatic rectal chylamdia.
-#' @param uct.dur.asympt Average in weeks duration of asymptomatic urethral chylamdia.
-#' @param ct.dur.tx Average in weeks duration of treated chylamdia (both sites).
-#' @param ct.dur.ntx Average in weeks duration of untreated, symptomatic chylamdia (both sites).
+#' @param rct.dur.asympt Average in weeks duration of asymptomatic rectal chlamydia.
+#' @param uct.dur.asympt Average in weeks duration of asymptomatic urethral chlamydia.
+#' @param ct.dur.tx Average in weeks duration of treated chlamydia (both sites).
+#' @param ct.dur.ntx Average in weeks duration of untreated, symptomatic chlamydia (both sites).
 #'        If \code{NULL}, uses site-specific durations for asymptomatic infections.
 #' @param gc.prob.cease Probability of ceasing sexual activity during symptomatic
 #'        infection with gonorrhea.
 #' @param ct.prob.cease Probability of ceasing sexual activity during symptomatic
-#'        infection with chylamdia.
+#'        infection with chlamydia.
 #' @param gc.sympt.prob.tx Probability of treatment for symptomatic gonorrhea.
-#' @param ct.sympt.prob.tx Probability of treatment for symptomatic chylamdia.
+#' @param ct.sympt.prob.tx Probability of treatment for symptomatic chlamydia.
 #' @param gc.asympt.prob.tx Probability of treatment for asymptomatic gonorrhea.
-#' @param ct.asympt.prob.tx Probability of treatment for asymptomatic chylamdia.
+#' @param ct.asympt.prob.tx Probability of treatment for asymptomatic chlamydia.
 #' @param prep.sti.screen.int Interval in days between STI screening at PrEP visits.
 #' @param prep.sti.prob.tx Probability of treatment given positive screening during
 #'        PrEP visit.
@@ -280,10 +280,10 @@
 #'        a condom used by the insertive partner.
 #' @param hiv.rgc.rr Relative risk of HIV infection given current rectal gonorrhea.
 #' @param hiv.ugc.rr Relative risk of HIV infection given current urethral gonorrhea.
-#' @param hiv.rct.rr Relative risk of HIV infection given current rectal chylamdia.
-#' @param hiv.uct.rr Relative risk of HIV infection given current urethral chylamdia.
+#' @param hiv.rct.rr Relative risk of HIV infection given current rectal chlamydia.
+#' @param hiv.uct.rr Relative risk of HIV infection given current urethral chlamydia.
 #' @param hiv.dual.rr Additive proportional risk, from 0 to 1, for HIV infection
-#'        given dual infection with both gonorrhea and chylamdia.
+#'        given dual infection with both gonorrhea and chlamydia.
 #'
 #' @param ... Additional arguments passed to the function.
 #'
@@ -555,6 +555,10 @@ param_msm <- function(nwstats,
 #'        \code{nwstats} output from \code{\link{calc_nwstats_msm}}.
 #' @param prev.B Initial disease prevalence among black MSM.
 #' @param prev.W Initial disease prevalence among white MSM.
+#' @param prev.ugc Initial prevalence of urethral gonorrhea.
+#' @param prev.rgc Initial prevalence of rectal gonorrhea.
+#' @param prev.uct Initial prevalence of urethral chlamydia.
+#' @param prev.rct Initial prevalence of rectal chlamydia.
 #' @param ... Additional arguments passed to function.
 #'
 #' @return
@@ -565,8 +569,12 @@ param_msm <- function(nwstats,
 #'
 #' @export
 init_msm <- function(nwstats,
-                     prev.B = 0.15,
-                     prev.W = 0.15,
+                     prev.B = 0.253,
+                     prev.W = 0.253,
+                     prev.ugc = 0.01,
+                     prev.rgc = 0.05,
+                     prev.uct = 0.01,
+                     prev.rct = 0.05,
                      ...) {
 
   p <- get_args(formal.args = formals(sys.function()),
diff --git a/man/init_msm.Rd b/man/init_msm.Rd
index e14f1d1c..8ba09a00 100644
--- a/man/init_msm.Rd
+++ b/man/init_msm.Rd
@@ -4,7 +4,8 @@
 \alias{init_msm}
 \title{Epidemic Model Initial Conditions}
 \usage{
-init_msm(nwstats, prev.B = 0.15, prev.W = 0.15, ...)
+init_msm(nwstats, prev.B = 0.253, prev.W = 0.253, prev.ugc = 0.01,
+  prev.rgc = 0.05, prev.uct = 0.01, prev.rct = 0.05, ...)
 }
 \arguments{
 \item{nwstats}{Target statistics for the network model. An object of class
@@ -14,6 +15,14 @@ init_msm(nwstats, prev.B = 0.15, prev.W = 0.15, ...)
 
 \item{prev.W}{Initial disease prevalence among white MSM.}
 
+\item{prev.ugc}{Initial prevalence of urethral gonorrhea.}
+
+\item{prev.rgc}{Initial prevalence of rectal gonorrhea.}
+
+\item{prev.uct}{Initial prevalence of urethral chlamydia.}
+
+\item{prev.rct}{Initial prevalence of rectal chlamydia.}
+
 \item{...}{Additional arguments passed to function.}
 }
 \value{
diff --git a/man/param_msm.Rd b/man/param_msm.Rd
index 39ac31eb..e6416968 100644
--- a/man/param_msm.Rd
+++ b/man/param_msm.Rd
@@ -406,9 +406,9 @@ partnerships only, versus all partnerships.}
 
 \item{ugc.tprob}{Probability of urethral gonorrhea infection per act.}
 
-\item{rct.tprob}{Probability of rectal chylamdia infection per act.}
+\item{rct.tprob}{Probability of rectal chlamydia infection per act.}
 
-\item{uct.tprob}{Probability of urethral chylamdia infection per act.}
+\item{uct.tprob}{Probability of urethral chlamydia infection per act.}
 
 \item{rgc.sympt.prob}{Probability of symptoms given infection with rectal
 gonorrhea.}
@@ -417,10 +417,10 @@ gonorrhea.}
 gonorrhea.}
 
 \item{rct.sympt.prob}{Probability of symptoms given infection with rectal
-chylamdia.}
+chlamydia.}
 
 \item{uct.sympt.prob}{Probability of symptoms given infection with urethral
-chylamdia.}
+chlamydia.}
 
 \item{rgc.dur.asympt}{Average duration in weeks of asymptomatic rectal gonorrhea.}
 
@@ -431,28 +431,28 @@ chylamdia.}
 \item{gc.dur.ntx}{Average duration in weeks of untreated, symptomatic gonorrhea (both sites).
 If \code{NULL}, uses site-specific durations for asymptomatic infections.}
 
-\item{rct.dur.asympt}{Average in weeks duration of asymptomatic rectal chylamdia.}
+\item{rct.dur.asympt}{Average in weeks duration of asymptomatic rectal chlamydia.}
 
-\item{uct.dur.asympt}{Average in weeks duration of asymptomatic urethral chylamdia.}
+\item{uct.dur.asympt}{Average in weeks duration of asymptomatic urethral chlamydia.}
 
-\item{ct.dur.tx}{Average in weeks duration of treated chylamdia (both sites).}
+\item{ct.dur.tx}{Average in weeks duration of treated chlamydia (both sites).}
 
-\item{ct.dur.ntx}{Average in weeks duration of untreated, symptomatic chylamdia (both sites).
+\item{ct.dur.ntx}{Average in weeks duration of untreated, symptomatic chlamydia (both sites).
 If \code{NULL}, uses site-specific durations for asymptomatic infections.}
 
 \item{gc.prob.cease}{Probability of ceasing sexual activity during symptomatic
 infection with gonorrhea.}
 
 \item{ct.prob.cease}{Probability of ceasing sexual activity during symptomatic
-infection with chylamdia.}
+infection with chlamydia.}
 
 \item{gc.sympt.prob.tx}{Probability of treatment for symptomatic gonorrhea.}
 
-\item{ct.sympt.prob.tx}{Probability of treatment for symptomatic chylamdia.}
+\item{ct.sympt.prob.tx}{Probability of treatment for symptomatic chlamydia.}
 
 \item{gc.asympt.prob.tx}{Probability of treatment for asymptomatic gonorrhea.}
 
-\item{ct.asympt.prob.tx}{Probability of treatment for asymptomatic chylamdia.}
+\item{ct.asympt.prob.tx}{Probability of treatment for asymptomatic chlamydia.}
 
 \item{prep.sti.screen.int}{Interval in days between STI screening at PrEP visits.}
 
@@ -469,12 +469,12 @@ a condom used by the insertive partner.}
 
 \item{hiv.ugc.rr}{Relative risk of HIV infection given current urethral gonorrhea.}
 
-\item{hiv.rct.rr}{Relative risk of HIV infection given current rectal chylamdia.}
+\item{hiv.rct.rr}{Relative risk of HIV infection given current rectal chlamydia.}
 
-\item{hiv.uct.rr}{Relative risk of HIV infection given current urethral chylamdia.}
+\item{hiv.uct.rr}{Relative risk of HIV infection given current urethral chlamydia.}
 
 \item{hiv.dual.rr}{Additive proportional risk, from 0 to 1, for HIV infection
-given dual infection with both gonorrhea and chylamdia.}
+given dual infection with both gonorrhea and chlamydia.}
 
 \item{...}{Additional arguments passed to the function.}
 }

From d255cefd186c712196afc81164cb0e956402894b Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Tue, 9 Aug 2016 10:07:58 -0400
Subject: [PATCH 03/40] Add EpiModelHPC to depends

---
 DESCRIPTION             | 1 +
 NAMESPACE               | 1 +
 R/EpiModelHIV-package.R | 2 +-
 3 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 37e1fa41..3bc43132 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -14,6 +14,7 @@ License: GPL-2
 Depends:
     R (>= 3.2.0),
     EpiModel (>= 1.2.7),
+    EpiModelHPC (>= 1.3.0),
     ergm (>= 3.5),
     tergm,
     tergmLite
diff --git a/NAMESPACE b/NAMESPACE
index d8e0a138..917b82da 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -63,6 +63,7 @@ export(verbose_msm)
 export(vl_het)
 export(vl_msm)
 import(EpiModel)
+import(EpiModelHPC)
 import(bindata)
 import(ergm)
 import(network)
diff --git a/R/EpiModelHIV-package.R b/R/EpiModelHIV-package.R
index 37b59fdb..edb7a5da 100644
--- a/R/EpiModelHIV-package.R
+++ b/R/EpiModelHIV-package.R
@@ -18,7 +18,7 @@
 #' @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

From c7d133beefab992000c398f941a91e2cb8bc4da7 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Tue, 9 Aug 2016 11:38:31 -0400
Subject: [PATCH 04/40] Move save.other internal

---
 R/EpiModelHIV-package.R    | 2 +-
 man/EpiModelHIV-package.Rd | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/EpiModelHIV-package.R b/R/EpiModelHIV-package.R
index edb7a5da..17373f76 100644
--- a/R/EpiModelHIV-package.R
+++ b/R/EpiModelHIV-package.R
@@ -5,7 +5,7 @@
 #'    Type: \tab Package\cr
 #'    Version: \tab 1.0.0\cr
 #'    Date: \tab 2016-06-25\cr
-#'    License: \tab GPL (>= 2)\cr
+#'    License: \tab GPL-3\cr
 #'    LazyLoad: \tab yes\cr
 #' }
 #'
diff --git a/man/EpiModelHIV-package.Rd b/man/EpiModelHIV-package.Rd
index 368bcba8..3daf0084 100644
--- a/man/EpiModelHIV-package.Rd
+++ b/man/EpiModelHIV-package.Rd
@@ -11,7 +11,7 @@
    Type: \tab Package\cr
    Version: \tab 1.0.0\cr
    Date: \tab 2016-06-25\cr
-   License: \tab GPL (>= 2)\cr
+   License: \tab GPL-3\cr
    LazyLoad: \tab yes\cr
 }
 }

From 62c5a1404ef7b7a0eecd6973ecc8ab78e82c21fb Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Tue, 9 Aug 2016 11:38:52 -0400
Subject: [PATCH 05/40] Update docs

---
 R/params.R         | 5 ++---
 man/control_msm.Rd | 7 ++-----
 2 files changed, 4 insertions(+), 8 deletions(-)

diff --git a/R/params.R b/R/params.R
index 0984e1d7..93c5ce79 100644
--- a/R/params.R
+++ b/R/params.R
@@ -643,8 +643,6 @@ init_msm <- function(nwstats,
 #'        external text files.
 #' @param save.nwstats Calculate and save network statistics as defined in the
 #'        \code{simnet} modules.
-#' @param save.other Character vector containing other list elements of \code{dat}
-#'        to save.
 #' @param verbose If \code{TRUE}, print out simulation progress to the console
 #'        if in interactive mode or text files if in batch mode.
 #' @param verbose.int Integer specifying the interval between time steps at which
@@ -687,7 +685,6 @@ control_msm <- function(simno = 1,
                         prev.FUN = prevalence_msm,
                         verbose.FUN = verbose_msm,
                         save.nwstats = FALSE,
-                        save.other = c("attr", "temp", "el", "p"),
                         verbose = TRUE,
                         verbose.int = 1,
                         ...) {
@@ -705,6 +702,8 @@ control_msm <- function(simno = 1,
   p$bi.mods <- bi.mods
   p$user.mods <- grep(".FUN", names(dot.args), value = TRUE)
 
+  p$save.other = c("attr", "temp", "el", "p")
+
   p$save.network = FALSE
 
   class(p) <- "control.net"
diff --git a/man/control_msm.Rd b/man/control_msm.Rd
index 3a033034..3b4ab862 100644
--- a/man/control_msm.Rd
+++ b/man/control_msm.Rd
@@ -14,8 +14,8 @@ control_msm(simno = 1, nsims = 1, ncores = 1, nsteps = 100, start = 1,
   riskhist.FUN = riskhist_msm, position.FUN = position_msm,
   trans.FUN = trans_msm, stitrans.FUN = sti_trans,
   stirecov.FUN = sti_recov, stitx.FUN = sti_tx, prev.FUN = prevalence_msm,
-  verbose.FUN = verbose_msm, save.nwstats = FALSE, save.other = c("attr",
-  "temp", "el", "p"), verbose = TRUE, verbose.int = 1, ...)
+  verbose.FUN = verbose_msm, save.nwstats = FALSE, verbose = TRUE,
+  verbose.int = 1, ...)
 }
 \arguments{
 \item{simno}{Unique ID for the simulation run, used for file naming purposes
@@ -90,9 +90,6 @@ external text files.}
 \item{save.nwstats}{Calculate and save network statistics as defined in the
 \code{simnet} modules.}
 
-\item{save.other}{Character vector containing other list elements of \code{dat}
-to save.}
-
 \item{verbose}{If \code{TRUE}, print out simulation progress to the console
 if in interactive mode or text files if in batch mode.}
 

From b45f0189e86e72f9b79054ba529a78bd27c8e156 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Tue, 9 Aug 2016 11:39:08 -0400
Subject: [PATCH 06/40] Bump HPC required version

---
 DESCRIPTION | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 3bc43132..3aab0fe7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -10,11 +10,11 @@ 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
   sub-Saharan Africa.
-License: GPL-2
+License: GPL-3
 Depends:
     R (>= 3.2.0),
     EpiModel (>= 1.2.7),
-    EpiModelHPC (>= 1.3.0),
+    EpiModelHPC (>= 1.3.1),
     ergm (>= 3.5),
     tergm,
     tergmLite

From 52dcdde9ae035d4795e1868ce3b75d8f396081d3 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Wed, 10 Aug 2016 10:34:38 -0400
Subject: [PATCH 07/40] Fix stop check in reinit_msm

---
 R/mod.initialize.R | 14 ++++++++++----
 man/reinit_msm.Rd  |  5 +++++
 2 files changed, 15 insertions(+), 4 deletions(-)

diff --git a/R/mod.initialize.R b/R/mod.initialize.R
index ce974408..7010216b 100644
--- a/R/mod.initialize.R
+++ b/R/mod.initialize.R
@@ -760,6 +760,11 @@ init_ccr5_msm <- function(dat) {
 #' @param x An \code{EpiModel} object of class \code{\link{netsim}}.
 #' @inheritParams initialize_msm
 #'
+#' @details
+#' Currently, the necessary components that must be on \code{x} for a simulation
+#' to be restarted must be: param, control, nwparam, epi, attr, temp, el, p.
+#' TODO: describe this more.
+#'
 #' @return
 #' This function resets the data elements on the \code{dat} master data object
 #' in the needed ways for the time loop to function.
@@ -769,10 +774,11 @@ init_ccr5_msm <- function(dat) {
 #'
 reinit_msm <- function(x, param, init, control, s) {
 
-  if (any(c("param", "control", "nwparam", "epi", "attr", "temp",
-    "el", "p") %in% names(x)) == FALSE) {
-    stop("x must contain the following elements for restarting: param control",
-         "nwparam epi attr temp el p", call. = FALSE)
+  need.for.reinit <- c("param", "control", "nwparam", "epi", "attr", "temp", "el", "p")
+  if (!all(need.for.reinit %in% names(x))) {
+    stop("x must contain the following elements for restarting: ",
+         "param, control, nwparam, epi, attr, temp, el, p",
+         call. = FALSE)
   }
 
   if (length(x$el) == 1) {
diff --git a/man/reinit_msm.Rd b/man/reinit_msm.Rd
index 312b3c82..1f526172 100644
--- a/man/reinit_msm.Rd
+++ b/man/reinit_msm.Rd
@@ -25,6 +25,11 @@ in the needed ways for the time loop to function.
 This function reinitializes an epidemic model to restart at a
              specified time step given an input \code{netsim} object.
 }
+\details{
+Currently, the necessary components that must be on \code{x} for a simulation
+to be restarted must be: param, control, nwparam, epi, attr, temp, el, p.
+TODO: describe this more.
+}
 \keyword{module}
 \keyword{msm}
 

From 6feb05674dca1db29d10e3f9ed9af1073c7c7196 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Wed, 10 Aug 2016 11:20:01 -0400
Subject: [PATCH 08/40] Use markdown for author table

---
 README.md | 20 +++++++-------------
 1 file changed, 7 insertions(+), 13 deletions(-)

diff --git a/README.md b/README.md
index 669a9e7d..8046f6ce 100644
--- a/README.md
+++ b/README.md
@@ -23,19 +23,13 @@ Documentation on using this software package is forthcoming, although limited fu
 ### Note on Repository and Package Name
 This Github repository `EpiModelHIV` and the R package within it were previously named `EpiModelHIVmsm`. On 6/24/2016, we merged that MSM package with our `EpiModelHIVhet` (HIV models for heterosexuals) into this combined repository and package. All scripts from those separate packages should still function with this current version after changing the input to `library`. 
 
-## Authors
-<table>
-  <tr>
-    <td><a href="http://samueljenness.org/" target="_blank">Samuel M. Jenness</a></th>
-    <td>Department of Epidemiology</th>
-    <td>Emory University</th>
-  </tr>
-  <tr>
-    <td><a href="http://faculty.washington.edu/goodreau/" target="_blank">Steven M. Goodreau</a></td>
-    <td>Department of Anthropology</td>
-    <td>University of Washington</td>
-  </tr>
-</table>
+## Software Development Team
+
+Author                                                        | Department                 | Institution
+------------------------------------------------------------- | -------------------------- | ------------------------
+[Samuel M. Jenness](http://samueljenness.org/)                | Department of Epidemiology | Emory University
+[Steven M. Goodreau](http://faculty.washington.edu/goodreau/) | Department of Anthropology | University of Washington
+
 
 ## Literature
 

From 59f4da59a81feb38046eeec5ead870108c9b6058 Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Wed, 10 Aug 2016 21:39:10 -0400
Subject: [PATCH 09/40] Cosmetic update

---
 R/mod.sti.R | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/R/mod.sti.R b/R/mod.sti.R
index 25fbcf9e..931acfb6 100644
--- a/R/mod.sti.R
+++ b/R/mod.sti.R
@@ -539,7 +539,8 @@ sti_tx <- function(dat, at) {
                                       dat$attr$rCT.infTime < at &
                                       is.na(dat$attr$rCT.tx)))
   idsUCT_prep_tx <- intersect(idsSTI_screen,
-                              which(dat$attr$uCT == 1 & dat$attr$uCT.infTime < at &
+                              which(dat$attr$uCT == 1 &
+                                      dat$attr$uCT.infTime < at &
                                       is.na(dat$attr$uCT.tx)))
 
   txRGC_prep <- idsRGC_prep_tx[which(rbinom(length(idsRGC_prep_tx), 1,

From 709621aa75e745c1c18d522a6d958b3c6792d063 Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Wed, 10 Aug 2016 21:39:57 -0400
Subject: [PATCH 10/40] Reset sti screening interval when stopping prep

---
 R/mod.prep.R | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/R/mod.prep.R b/R/mod.prep.R
index 44afed97..fb97b5c6 100644
--- a/R/mod.prep.R
+++ b/R/mod.prep.R
@@ -28,6 +28,7 @@ prep_msm <- function(dat, at) {
   prepClass <- dat$attr$prepClass
   prepLastRisk <- dat$attr$prepLastRisk
   prepStartTime <- dat$attr$prepStartTime
+  prepLastStiScreen <- dat$attr$prepLastStiScreen
 
   # Parameters
 
@@ -78,6 +79,7 @@ prep_msm <- function(dat, at) {
   prepStat[idsStp] <- 0
   prepLastRisk[idsStp] <- NA
   prepStartTime[idsStp] <- NA
+  prepLastStiScreen[idsStp] <- NA
 
 
   ## Initiation ----------------------------------------------------------------
@@ -120,6 +122,7 @@ prep_msm <- function(dat, at) {
   dat$attr$prepStartTime <- prepStartTime
   dat$attr$prepClass <- prepClass
   dat$attr$prepLastRisk <- prepLastRisk
+  dat$attr$prepLastStiScreen <- prepLastStiScreen
 
   # Summary Statistics
   dat$epi$prepCov[at] <- prepCov

From 111c3aac39c742cacf7fe3baacc059b550d8ad13 Mon Sep 17 00:00:00 2001
From: kweiss2 <kevin.weiss@emory.edu>
Date: Mon, 15 Aug 2016 12:01:09 -0400
Subject: [PATCH 11/40] Combined STI prev

---
 R/mod.prevalence.R | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/R/mod.prevalence.R b/R/mod.prevalence.R
index 3ea3e8fe..ca5e923c 100644
--- a/R/mod.prevalence.R
+++ b/R/mod.prevalence.R
@@ -156,6 +156,10 @@ prevalence_msm <- function(dat, at) {
   dat$epi$ir100.uct[at] <- (dat$epi$incid.uct[at] / sum(uCT == 0, na.rm = TRUE)) * 5200
   dat$epi$ir100.ct[at] <- (dat$epi$incid.ct[at]/ sum(rCT == 0 | uCT == 0, na.rm = TRUE)) * 5200
 
+  dat$epi$prev.sti[at] <- sum(rGC == 1 | ugc == 1 | rCT ==1 | uCT == 1, na.rm = TRUE) / dat$epi$num[at]
+  dat$epi$ir100.sti[at] <- ((dat$epi$incid.ct[at] + dat$epi$incid.gc)/ sum(rCT == 0 | uCT == 0 | rGC == 0 | uGC ==0, na.rm = TRUE)) * 5200
+  
+  
   return(dat)
 }
 

From d5ed91fd9ac8658cada35f36fead11a619a31be1 Mon Sep 17 00:00:00 2001
From: kweiss2 <kevin.weiss@emory.edu>
Date: Mon, 15 Aug 2016 13:16:00 -0400
Subject: [PATCH 12/40] Typo fix

---
 R/mod.prevalence.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/mod.prevalence.R b/R/mod.prevalence.R
index ca5e923c..6e307c0b 100644
--- a/R/mod.prevalence.R
+++ b/R/mod.prevalence.R
@@ -156,8 +156,8 @@ prevalence_msm <- function(dat, at) {
   dat$epi$ir100.uct[at] <- (dat$epi$incid.uct[at] / sum(uCT == 0, na.rm = TRUE)) * 5200
   dat$epi$ir100.ct[at] <- (dat$epi$incid.ct[at]/ sum(rCT == 0 | uCT == 0, na.rm = TRUE)) * 5200
 
-  dat$epi$prev.sti[at] <- sum(rGC == 1 | ugc == 1 | rCT ==1 | uCT == 1, na.rm = TRUE) / dat$epi$num[at]
-  dat$epi$ir100.sti[at] <- ((dat$epi$incid.ct[at] + dat$epi$incid.gc)/ sum(rCT == 0 | uCT == 0 | rGC == 0 | uGC ==0, na.rm = TRUE)) * 5200
+  dat$epi$prev.sti[at] <- sum(rGC == 1 | uGC == 1 | rCT ==1 | uCT == 1, na.rm = TRUE) / dat$epi$num[at]
+  dat$epi$ir100.sti[at] <- ((dat$epi$incid.ct[at] + dat$epi$incid.gc)/ sum(rCT == 0 | uCT == 0 | rGC == 0 | uGC == 0, na.rm = TRUE)) * 5200
   
   
   return(dat)

From 6da9e55febda515e5f6857908ada78f545d79327 Mon Sep 17 00:00:00 2001
From: kweiss2 <kevin.weiss@emory.edu>
Date: Mon, 15 Aug 2016 13:20:51 -0400
Subject: [PATCH 13/40] Fix to STI incidence statistics

---
 R/mod.prevalence.R | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/R/mod.prevalence.R b/R/mod.prevalence.R
index 6e307c0b..b7e006ae 100644
--- a/R/mod.prevalence.R
+++ b/R/mod.prevalence.R
@@ -150,14 +150,14 @@ prevalence_msm <- function(dat, at) {
 
   dat$epi$ir100.rgc[at] <- (dat$epi$incid.rgc[at] / sum(rGC == 0, na.rm = TRUE)) * 5200
   dat$epi$ir100.ugc[at] <- (dat$epi$incid.ugc[at] / sum(uGC == 0, na.rm = TRUE)) * 5200
-  dat$epi$ir100.gc[at] <- (dat$epi$incid.gc[at]/ sum(rGC == 0 | uGC == 0, na.rm = TRUE)) * 5200
+  dat$epi$ir100.gc[at] <- (dat$epi$incid.gc[at]/ sum(rGC == 0 & uGC == 0, na.rm = TRUE)) * 5200
 
   dat$epi$ir100.rct[at] <- (dat$epi$incid.rct[at] / sum(rCT == 0, na.rm = TRUE)) * 5200
   dat$epi$ir100.uct[at] <- (dat$epi$incid.uct[at] / sum(uCT == 0, na.rm = TRUE)) * 5200
-  dat$epi$ir100.ct[at] <- (dat$epi$incid.ct[at]/ sum(rCT == 0 | uCT == 0, na.rm = TRUE)) * 5200
+  dat$epi$ir100.ct[at] <- (dat$epi$incid.ct[at]/ sum(rCT == 0 & uCT == 0, na.rm = TRUE)) * 5200
 
   dat$epi$prev.sti[at] <- sum(rGC == 1 | uGC == 1 | rCT ==1 | uCT == 1, na.rm = TRUE) / dat$epi$num[at]
-  dat$epi$ir100.sti[at] <- ((dat$epi$incid.ct[at] + dat$epi$incid.gc)/ sum(rCT == 0 | uCT == 0 | rGC == 0 | uGC == 0, na.rm = TRUE)) * 5200
+  dat$epi$ir100.sti[at] <- ((dat$epi$incid.ct[at] + dat$epi$incid.gc)/ sum(rCT == 0 & uCT == 0 & rGC == 0 & uGC == 0, na.rm = TRUE)) * 5200
   
   
   return(dat)

From a718cbf4a9a7c85d783dcd2f61dc4191ac714756 Mon Sep 17 00:00:00 2001
From: kweiss2 <kevin.weiss@emory.edu>
Date: Mon, 15 Aug 2016 15:48:48 -0400
Subject: [PATCH 14/40] Attempted fix

---
 R/mod.prevalence.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/mod.prevalence.R b/R/mod.prevalence.R
index b7e006ae..8a022215 100644
--- a/R/mod.prevalence.R
+++ b/R/mod.prevalence.R
@@ -157,7 +157,7 @@ prevalence_msm <- function(dat, at) {
   dat$epi$ir100.ct[at] <- (dat$epi$incid.ct[at]/ sum(rCT == 0 & uCT == 0, na.rm = TRUE)) * 5200
 
   dat$epi$prev.sti[at] <- sum(rGC == 1 | uGC == 1 | rCT ==1 | uCT == 1, na.rm = TRUE) / dat$epi$num[at]
-  dat$epi$ir100.sti[at] <- ((dat$epi$incid.ct[at] + dat$epi$incid.gc)/ sum(rCT == 0 & uCT == 0 & rGC == 0 & uGC == 0, na.rm = TRUE)) * 5200
+  dat$epi$ir100.sti[at] <- ((dat$epi$incid.ct[at] + dat$epi$incid.gc[at])/ sum(rCT == 0 & uCT == 0 & rGC == 0 & uGC == 0, na.rm = TRUE)) * 5200
   
   
   return(dat)

From 206da51ddaf8bb639be2989b2cd89fb98c1458ca Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Wed, 17 Aug 2016 22:43:16 -0400
Subject: [PATCH 15/40] Fix error in calculating combined STI incidence

---
 R/mod.prevalence.R | 21 +++++++++++++++------
 1 file changed, 15 insertions(+), 6 deletions(-)

diff --git a/R/mod.prevalence.R b/R/mod.prevalence.R
index 8a022215..51d02292 100644
--- a/R/mod.prevalence.R
+++ b/R/mod.prevalence.R
@@ -150,16 +150,25 @@ prevalence_msm <- function(dat, at) {
 
   dat$epi$ir100.rgc[at] <- (dat$epi$incid.rgc[at] / sum(rGC == 0, na.rm = TRUE)) * 5200
   dat$epi$ir100.ugc[at] <- (dat$epi$incid.ugc[at] / sum(uGC == 0, na.rm = TRUE)) * 5200
-  dat$epi$ir100.gc[at] <- (dat$epi$incid.gc[at]/ sum(rGC == 0 & uGC == 0, na.rm = TRUE)) * 5200
+  dat$epi$ir100.gc[at] <- (dat$epi$incid.gc[at] /
+                             (sum(rGC == 0, na.rm = TRUE) +
+                                sum(uGC == 0, na.rm = TRUE))) * 5200
 
   dat$epi$ir100.rct[at] <- (dat$epi$incid.rct[at] / sum(rCT == 0, na.rm = TRUE)) * 5200
   dat$epi$ir100.uct[at] <- (dat$epi$incid.uct[at] / sum(uCT == 0, na.rm = TRUE)) * 5200
-  dat$epi$ir100.ct[at] <- (dat$epi$incid.ct[at]/ sum(rCT == 0 & uCT == 0, na.rm = TRUE)) * 5200
+  dat$epi$ir100.ct[at] <- (dat$epi$incid.ct[at] /
+                             (sum(rCT == 0, na.rm = TRUE) +
+                                sum(uCT == 0, na.rm = TRUE))) * 5200
+
+  dat$epi$prev.sti[at] <- sum(rGC == 1 | uGC == 1 |
+                                rCT ==1 | uCT == 1, na.rm = TRUE) / dat$epi$num[at]
+  dat$epi$ir100.sti[at] <- ((dat$epi$incid.ct[at] + dat$epi$incid.gc[at]) /
+                              (sum(rGC == 0, na.rm = TRUE) +
+                                 sum(uGC == 0, na.rm = TRUE) +
+                                 sum(rCT == 0, na.rm = TRUE) +
+                                 sum(uCT == 0, na.rm = TRUE))) * 5200
+
 
-  dat$epi$prev.sti[at] <- sum(rGC == 1 | uGC == 1 | rCT ==1 | uCT == 1, na.rm = TRUE) / dat$epi$num[at]
-  dat$epi$ir100.sti[at] <- ((dat$epi$incid.ct[at] + dat$epi$incid.gc[at])/ sum(rCT == 0 & uCT == 0 & rGC == 0 & uGC == 0, na.rm = TRUE)) * 5200
-  
-  
   return(dat)
 }
 

From f29579b5869dd44c814f782fbaab12db58cf150d Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Wed, 7 Sep 2016 23:47:32 -0400
Subject: [PATCH 16/40] Rcpp package related updates

---
 R/RcppExports.R     |  2 +-
 src/RcppExports.cpp | 18 +++++++++---------
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/R/RcppExports.R b/R/RcppExports.R
index d3682776..29b12d23 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -1,4 +1,4 @@
-# This file was generated by Rcpp::compileAttributes
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 #' @title Aging Module
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 7d2ad9b6..b5ca1d1a 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -1,4 +1,4 @@
-// This file was generated by Rcpp::compileAttributes
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
 // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 #include <Rcpp.h>
@@ -9,23 +9,23 @@ using namespace Rcpp;
 List aging_msm(List dat, int at);
 RcppExport SEXP EpiModelHIV_aging_msm(SEXP datSEXP, SEXP atSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< List >::type dat(datSEXP);
     Rcpp::traits::input_parameter< int >::type at(atSEXP);
-    __result = Rcpp::wrap(aging_msm(dat, at));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(aging_msm(dat, at));
+    return rcpp_result_gen;
 END_RCPP
 }
 // aging_het
 List aging_het(List dat, int at);
 RcppExport SEXP EpiModelHIV_aging_het(SEXP datSEXP, SEXP atSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< List >::type dat(datSEXP);
     Rcpp::traits::input_parameter< int >::type at(atSEXP);
-    __result = Rcpp::wrap(aging_het(dat, at));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(aging_het(dat, at));
+    return rcpp_result_gen;
 END_RCPP
 }

From b999a476f1273c3247f25448f41ae2f9aee0ea60 Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Wed, 7 Sep 2016 23:48:01 -0400
Subject: [PATCH 17/40] Fix bug for dual-site infection

---
 R/mod.sti.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/mod.sti.R b/R/mod.sti.R
index 931acfb6..2ae8462b 100644
--- a/R/mod.sti.R
+++ b/R/mod.sti.R
@@ -573,10 +573,10 @@ sti_tx <- function(dat, at) {
 
   # add tx at other site
   dat$attr$rGC.tx[which(dat$attr$uGC.tx == 1 & dat$attr$rGC == 1)] <- 1
-  dat$attr$uGC.tx[which(dat$attr$rGC.tx == 1 & dat$attr$rGC == 1)] <- 1
+  dat$attr$uGC.tx[which(dat$attr$rGC.tx == 1 & dat$attr$uGC == 1)] <- 1
 
   dat$attr$rCT.tx[which(dat$attr$uCT.tx == 1 & dat$attr$rCT == 1)] <- 1
-  dat$attr$uCT.tx[which(dat$attr$rCT.tx == 1 & dat$attr$rCT == 1)] <- 1
+  dat$attr$uCT.tx[which(dat$attr$rCT.tx == 1 & dat$attr$uCT == 1)] <- 1
 
   # summary stats
   if (is.null(dat$epi$txGC)) {

From 22ff96a2df8c13f506ba095655e72185b013d697 Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Wed, 7 Sep 2016 23:49:02 -0400
Subject: [PATCH 18/40] Fix bugs in STI recovery module

---
 R/mod.sti.R | 138 ++++++++++++++++++++++++++++++----------------------
 1 file changed, 79 insertions(+), 59 deletions(-)

diff --git a/R/mod.sti.R b/R/mod.sti.R
index 2ae8462b..7f5e3ab0 100644
--- a/R/mod.sti.R
+++ b/R/mod.sti.R
@@ -309,43 +309,53 @@ sti_recov <- function(dat, at) {
 
 
   # GC recovery
-  idsRGC_asympt <- which(dat$attr$rGC == 1 & dat$attr$rGC.infTime < at &
-                           dat$attr$rGC.sympt == 0)
-  idsUGC_asympt <- which(dat$attr$uGC == 1 & dat$attr$uGC.infTime < at &
-                           dat$attr$uGC.sympt == 0)
-  idsRGC_tx <- which(dat$attr$rGC == 1 & dat$attr$rGC.infTime < at &
-                       dat$attr$rGC.sympt == 1 & dat$attr$rGC.tx == 1)
-  idsUGC_tx <- which(dat$attr$uGC == 1 & dat$attr$uGC.infTime < at &
-                       dat$attr$uGC.sympt == 1 & dat$attr$uGC.tx == 1)
-  idsRGC_ntx <- which(dat$attr$rGC == 1 & dat$attr$rGC.infTime < at &
-                        dat$attr$rGC.sympt == 1 & dat$attr$rGC.tx == 0)
-  idsUGC_ntx <- which(dat$attr$uGC == 1 & dat$attr$uGC.infTime < at &
-                        dat$attr$uGC.sympt == 1 & dat$attr$uGC.tx == 0)
-
-  recovRGC_asympt <- idsRGC_asympt[which(rbinom(length(idsRGC_asympt), 1,
-                                                1/rgc.dur.asympt) == 1)]
-  recovUGC_asympt <- idsUGC_asympt[which(rbinom(length(idsUGC_asympt), 1,
-                                                1/ugc.dur.asympt) == 1)]
+  idsRGC_asympt_ntx <- which(dat$attr$rGC == 1 &
+                             dat$attr$rGC.infTime < at &
+                             dat$attr$rGC.sympt == 0 &
+                             (is.na(dat$attr$rGC.tx) | dat$attr$rGC.tx == 0))
+  idsUGC_asympt_ntx <- which(dat$attr$uGC == 1 &
+                             dat$attr$uGC.infTime < at &
+                             dat$attr$uGC.sympt == 0 &
+                             (is.na(dat$attr$uGC.tx) | dat$attr$uGC.tx == 0))
+  idsRGC_sympt_ntx <- which(dat$attr$rGC == 1 &
+                            dat$attr$rGC.infTime < at &
+                            dat$attr$rGC.sympt == 1 &
+                            (is.na(dat$attr$rGC.tx) | dat$attr$rGC.tx == 0))
+  idsUGC_sympt_ntx <- which(dat$attr$uGC == 1 &
+                            dat$attr$uGC.infTime < at &
+                            dat$attr$uGC.sympt == 1 &
+                            (is.na(dat$attr$uGC.tx) | dat$attr$uGC.tx == 0))
+  idsRGC_tx <- which(dat$attr$rGC == 1 &
+                     dat$attr$rGC.infTime < at &
+                     dat$attr$rGC.tx == 1)
+  idsUGC_tx <- which(dat$attr$uGC == 1 &
+                     dat$attr$uGC.infTime < at &
+                     dat$attr$uGC.tx == 1)
+
+  recovRGC_asympt_ntx <- idsRGC_asympt_ntx[which(rbinom(length(idsRGC_asympt_ntx), 1,
+                                                 1/rgc.dur.asympt) == 1)]
+  recovUGC_asympt_ntx <- idsUGC_asympt_ntx[which(rbinom(length(idsUGC_asympt_ntx), 1,
+                                                 1/ugc.dur.asympt) == 1)]
+
+  if (!is.null(gc.dur.ntx)) {
+    recovRGC_sympt_ntx <- idsRGC_sympt_ntx[which(rbinom(length(idsRGC_sympt_ntx), 1,
+                                                 1/gc.dur.ntx) == 1)]
+    recovUGC_sympt_ntx <- idsUGC_sympt_ntx[which(rbinom(length(idsUGC_sympt_ntx), 1,
+                                                 1/gc.dur.ntx) == 1)]
+  } else {
+    recovRGC_sympt_ntx <- idsRGC_sympt_ntx[which(rbinom(length(idsRGC_sympt_ntx), 1,
+                                                 1/rgc.dur.asympt) == 1)]
+    recovUGC_sympt_ntx <- idsUGC_sympt_ntx[which(rbinom(length(idsUGC_sympt_ntx), 1,
+                                                 1/ugc.dur.asympt) == 1)]
+  }
 
   recovRGC_tx <- idsRGC_tx[which(rbinom(length(idsRGC_tx), 1,
                                         1/gc.dur.tx) == 1)]
   recovUGC_tx <- idsUGC_tx[which(rbinom(length(idsUGC_tx), 1,
                                         1/gc.dur.tx) == 1)]
 
-  if (!is.null(gc.dur.ntx)) {
-    recovRGC_ntx <- idsRGC_ntx[which(rbinom(length(idsRGC_ntx), 1,
-                                            1/gc.dur.ntx) == 1)]
-    recovUGC_ntx <- idsUGC_ntx[which(rbinom(length(idsUGC_ntx), 1,
-                                            1/gc.dur.ntx) == 1)]
-  } else {
-    recovRGC_ntx <- idsRGC_ntx[which(rbinom(length(idsRGC_ntx), 1,
-                                            1/rgc.dur.asympt) == 1)]
-    recovUGC_ntx <- idsUGC_ntx[which(rbinom(length(idsUGC_ntx), 1,
-                                            1/ugc.dur.asympt) == 1)]
-  }
-
-  recovRGC <- c(recovRGC_asympt, recovRGC_tx, recovRGC_ntx)
-  recovUGC <- c(recovUGC_asympt, recovUGC_tx, recovUGC_ntx)
+  recovRGC <- c(recovRGC_asympt_ntx, recovRGC_sympt_ntx, recovRGC_tx)
+  recovUGC <- c(recovUGC_asympt_ntx, recovUGC_sympt_ntx, recovUGC_tx)
 
   dat$attr$rGC[recovRGC] <- 0
   dat$attr$rGC.sympt[recovRGC] <- NA
@@ -360,43 +370,53 @@ sti_recov <- function(dat, at) {
   dat$attr$GC.cease[c(recovRGC, recovUGC)] <- NA
 
   # CT recovery
-  idsRCT_asympt <- which(dat$attr$rCT == 1 & dat$attr$rCT.infTime < at &
-                           dat$attr$rCT.sympt == 0)
-  idsUCT_asympt <- which(dat$attr$uCT == 1 & dat$attr$uCT.infTime < at &
-                           dat$attr$uCT.sympt == 0)
-  idsRCT_tx <- which(dat$attr$rCT == 1 & dat$attr$rCT.infTime < at &
-                       dat$attr$rCT.sympt == 1 & dat$attr$rCT.tx == 1)
-  idsUCT_tx <- which(dat$attr$uCT == 1 & dat$attr$uCT.infTime < at &
-                       dat$attr$uCT.sympt == 1 & dat$attr$uCT.tx == 1)
-  idsRCT_ntx <- which(dat$attr$rCT == 1 & dat$attr$rCT.infTime < at &
-                        dat$attr$rCT.sympt == 1 & dat$attr$rCT.tx == 0)
-  idsUCT_ntx <- which(dat$attr$uCT == 1 & dat$attr$uCT.infTime < at &
-                        dat$attr$uCT.sympt == 1 & dat$attr$uCT.tx == 0)
-
-  recovRCT_asympt <- idsRCT_asympt[which(rbinom(length(idsRCT_asympt),
-                                                1, 1/rct.dur.asympt) == 1)]
-  recovUCT_asympt <- idsUCT_asympt[which(rbinom(length(idsUCT_asympt),
-                                                1, 1/uct.dur.asympt) == 1)]
-
-  recovRCT_tx <- idsRCT_tx[which(rbinom(length(idsRCT_tx),
-                                        1, 1/ct.dur.tx) == 1)]
-  recovUCT_tx <- idsUCT_tx[which(rbinom(length(idsUCT_tx),
-                                        1, 1/ct.dur.tx) == 1)]
+  idsRCT_asympt_ntx <- which(dat$attr$rCT == 1 &
+                             dat$attr$rCT.infTime < at &
+                             dat$attr$rCT.sympt == 0 &
+                             (is.na(dat$attr$rCT.tx) | dat$attr$rCT.tx == 0))
+  idsUCT_asympt_ntx <- which(dat$attr$uCT == 1 &
+                             dat$attr$uCT.infTime < at &
+                             dat$attr$uCT.sympt == 0 &
+                             (is.na(dat$attr$uCT.tx) | dat$attr$uCT.tx == 0))
+  idsRCT_sympt_ntx <- which(dat$attr$rCT == 1 &
+                            dat$attr$rCT.infTime < at &
+                            dat$attr$rCT.sympt == 1 &
+                            (is.na(dat$attr$rCT.tx) | dat$attr$rCT.tx == 0))
+  idsUCT_sympt_ntx <- which(dat$attr$uCT == 1 &
+                            dat$attr$uCT.infTime < at &
+                            dat$attr$uCT.sympt == 1 &
+                            (is.na(dat$attr$uCT.tx) | dat$attr$uCT.tx == 0))
+  idsRCT_tx <- which(dat$attr$rCT == 1 &
+                     dat$attr$rCT.infTime < at &
+                     dat$attr$rCT.tx == 1)
+  idsUCT_tx <- which(dat$attr$uCT == 1 &
+                     dat$attr$uCT.infTime < at &
+                     dat$attr$uCT.tx == 1)
+
+  recovRCT_asympt_ntx <- idsRCT_asympt_ntx[which(rbinom(length(idsRCT_asympt_ntx),
+                                                 1, 1/rct.dur.asympt) == 1)]
+  recovUCT_asympt_ntx <- idsUCT_asympt_ntx[which(rbinom(length(idsUCT_asympt_ntx),
+                                                 1, 1/uct.dur.asympt) == 1)]
 
   if (!is.null(ct.dur.ntx)) {
-    recovRCT_ntx <- idsRCT_ntx[which(rbinom(length(idsRCT_ntx),
+    recovRCT_sympt_ntx <- idsRCT_sympt_ntx[which(rbinom(length(idsRCT_sympt_ntx),
                                             1, 1/ct.dur.ntx) == 1)]
-    recovUCT_ntx <- idsUCT_ntx[which(rbinom(length(idsUCT_ntx),
+    recovUCT_sympt_ntx <- idsUCT_sympt_ntx[which(rbinom(length(idsUCT_sympt_ntx),
                                             1, 1/ct.dur.ntx) == 1)]
   } else {
-    recovRCT_ntx <- idsRCT_ntx[which(rbinom(length(idsRCT_ntx),
+    recovRCT_sympt_ntx <- idsRCT_sympt_ntx[which(rbinom(length(idsRCT_sympt_ntx),
                                             1, 1/rct.dur.asympt) == 1)]
-    recovUCT_ntx <- idsUCT_ntx[which(rbinom(length(idsUCT_ntx),
+    recovUCT_sympt_ntx <- idsUCT_sympt_ntx[which(rbinom(length(idsUCT_sympt_ntx),
                                             1, 1/uct.dur.asympt) == 1)]
   }
 
-  recovRCT <- c(recovRCT_asympt, recovRCT_tx, recovRCT_ntx)
-  recovUCT <- c(recovUCT_asympt, recovUCT_tx, recovUCT_ntx)
+  recovRCT_tx <- idsRCT_tx[which(rbinom(length(idsRCT_tx),
+                                        1, 1/ct.dur.tx) == 1)]
+  recovUCT_tx <- idsUCT_tx[which(rbinom(length(idsUCT_tx),
+                                        1, 1/ct.dur.tx) == 1)]
+
+  recovRCT <- c(recovRCT_asympt_ntx, recovRCT_sympt_ntx, recovRCT_tx)
+  recovUCT <- c(recovUCT_asympt_ntx, recovUCT_sympt_ntx, recovUCT_tx)
 
   dat$attr$rCT[recovRCT] <- 0
   dat$attr$rCT.sympt[recovRCT] <- NA

From d3650974401d91ad5bc784204b6467e81ec8dc08 Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Wed, 7 Sep 2016 23:49:24 -0400
Subject: [PATCH 19/40] Fix bugs in the STI tx module

---
 R/mod.sti.R | 74 +++++++++++++++++++++++++++++++++--------------------
 1 file changed, 46 insertions(+), 28 deletions(-)

diff --git a/R/mod.sti.R b/R/mod.sti.R
index 7f5e3ab0..e05d769d 100644
--- a/R/mod.sti.R
+++ b/R/mod.sti.R
@@ -470,12 +470,16 @@ sti_tx <- function(dat, at) {
   }
 
   # symptomatic gc treatment
-  idsRGC_tx_sympt <- which(dat$attr$rGC == 1 & dat$attr$rGC.infTime < at &
-                             dat$attr$rGC.sympt == 1 & is.na(dat$attr$rGC.tx) &
-                             dat$attr$prepStat %in% prep.stand.tx.grp)
-  idsUGC_tx_sympt <- which(dat$attr$uGC == 1 & dat$attr$uGC.infTime < at &
-                             dat$attr$uGC.sympt == 1 & is.na(dat$attr$uGC.tx) &
-                             dat$attr$prepStat %in% prep.stand.tx.grp)
+  idsRGC_tx_sympt <- which(dat$attr$rGC == 1 &
+                           dat$attr$rGC.infTime < at &
+                           dat$attr$rGC.sympt == 1 &
+                           is.na(dat$attr$rGC.tx) &
+                           dat$attr$prepStat %in% prep.stand.tx.grp)
+  idsUGC_tx_sympt <- which(dat$attr$uGC == 1 &
+                           dat$attr$uGC.infTime < at &
+                           dat$attr$uGC.sympt == 1 &
+                           is.na(dat$attr$uGC.tx) &
+                           dat$attr$prepStat %in% prep.stand.tx.grp)
   idsGC_tx_sympt <- c(idsRGC_tx_sympt, idsUGC_tx_sympt)
 
   txGC_sympt <- idsGC_tx_sympt[which(rbinom(length(idsGC_tx_sympt), 1,
@@ -484,12 +488,16 @@ sti_tx <- function(dat, at) {
   txUGC_sympt <- intersect(idsUGC_tx_sympt, txGC_sympt)
 
   # asymptomatic gc treatment
-  idsRGC_tx_asympt <- which(dat$attr$rGC == 1 & dat$attr$rGC.infTime < at &
-                              dat$attr$rGC.sympt == 0 & is.na(dat$attr$rGC.tx) &
-                              dat$attr$prepStat %in% prep.stand.tx.grp)
-  idsUGC_tx_asympt <- which(dat$attr$uGC == 1 & dat$attr$uGC.infTime < at &
-                              dat$attr$uGC.sympt == 0 & is.na(dat$attr$uGC.tx) &
-                              dat$attr$prepStat %in% prep.stand.tx.grp)
+  idsRGC_tx_asympt <- which(dat$attr$rGC == 1 &
+                            dat$attr$rGC.infTime < at &
+                            dat$attr$rGC.sympt == 0 &
+                            is.na(dat$attr$rGC.tx) &
+                            dat$attr$prepStat %in% prep.stand.tx.grp)
+  idsUGC_tx_asympt <- which(dat$attr$uGC == 1 &
+                            dat$attr$uGC.infTime < at &
+                            dat$attr$uGC.sympt == 0 &
+                            is.na(dat$attr$uGC.tx) &
+                            dat$attr$prepStat %in% prep.stand.tx.grp)
   idsGC_tx_asympt <- c(idsRGC_tx_asympt, idsUGC_tx_asympt)
 
   txGC_asympt <- idsGC_tx_asympt[which(rbinom(length(idsGC_tx_asympt), 1,
@@ -506,12 +514,16 @@ sti_tx <- function(dat, at) {
 
 
   # symptomatic ct treatment
-  idsRCT_tx_sympt <- which(dat$attr$rCT == 1 & dat$attr$rCT.infTime < at &
-                             dat$attr$rCT.sympt == 1 & is.na(dat$attr$rCT.tx) &
-                             dat$attr$prepStat %in% prep.stand.tx.grp)
-  idsUCT_tx_sympt <- which(dat$attr$uCT == 1 & dat$attr$uCT.infTime < at &
-                             dat$attr$uCT.sympt == 1 & is.na(dat$attr$uCT.tx) &
-                             dat$attr$prepStat %in% prep.stand.tx.grp)
+  idsRCT_tx_sympt <- which(dat$attr$rCT == 1 &
+                           dat$attr$rCT.infTime < at &
+                           dat$attr$rCT.sympt == 1 &
+                           is.na(dat$attr$rCT.tx) &
+                           dat$attr$prepStat %in% prep.stand.tx.grp)
+  idsUCT_tx_sympt <- which(dat$attr$uCT == 1 &
+                           dat$attr$uCT.infTime < at &
+                           dat$attr$uCT.sympt == 1 &
+                           is.na(dat$attr$uCT.tx) &
+                           dat$attr$prepStat %in% prep.stand.tx.grp)
   idsCT_tx_sympt <- c(idsRCT_tx_sympt, idsUCT_tx_sympt)
 
   txCT_sympt <- idsCT_tx_sympt[which(rbinom(length(idsCT_tx_sympt), 1,
@@ -520,12 +532,16 @@ sti_tx <- function(dat, at) {
   txUCT_sympt <- intersect(idsUCT_tx_sympt, txCT_sympt)
 
   # asymptomatic ct treatment
-  idsRCT_tx_asympt <- which(dat$attr$rCT == 1 & dat$attr$rCT.infTime < at &
-                              dat$attr$rCT.sympt == 0 & is.na(dat$attr$rCT.tx) &
-                              dat$attr$prepStat == 0)
-  idsUCT_tx_asympt <- which(dat$attr$uCT == 1 & dat$attr$uCT.infTime < at &
-                              dat$attr$uCT.sympt == 0 & is.na(dat$attr$uCT.tx) &
-                              dat$attr$prepStat == 0)
+  idsRCT_tx_asympt <- which(dat$attr$rCT == 1 &
+                            dat$attr$rCT.infTime < at &
+                            dat$attr$rCT.sympt == 0 &
+                            is.na(dat$attr$rCT.tx) &
+                            dat$attr$prepStat == 0)
+  idsUCT_tx_asympt <- which(dat$attr$uCT == 1 &
+                            dat$attr$uCT.infTime < at &
+                            dat$attr$uCT.sympt == 0 &
+                            is.na(dat$attr$uCT.tx) &
+                            dat$attr$prepStat == 0)
   idsCT_tx_asympt <- c(idsRCT_tx_asympt, idsUCT_tx_asympt)
 
   txCT_asympt <- idsCT_tx_asympt[which(rbinom(length(idsCT_tx_asympt), 1,
@@ -546,22 +562,24 @@ sti_tx <- function(dat, at) {
 
   dat$attr$prepLastStiScreen[idsSTI_screen] <- at
 
+  ## TODO: current approach does not allow for fixed prep-specific non-treatment
+  ##       attribute, where the prep.sti.prob.tx parameter < 1 (default)
   idsRGC_prep_tx <- intersect(idsSTI_screen,
                               which(dat$attr$rGC == 1 &
                                       dat$attr$rGC.infTime < at &
-                                      is.na(dat$attr$rGC.tx)))
+                                      (is.na(dat$attr$rGC.tx) | dat$attr$rGC.tx == 0)))
   idsUGC_prep_tx <- intersect(idsSTI_screen,
                               which(dat$attr$uGC == 1 &
                                       dat$attr$uGC.infTime < at &
-                                      is.na(dat$attr$uGC.tx)))
+                                      (is.na(dat$attr$uGC.tx) | dat$attr$uGC.tx == 0)))
   idsRCT_prep_tx <- intersect(idsSTI_screen,
                               which(dat$attr$rCT == 1 &
                                       dat$attr$rCT.infTime < at &
-                                      is.na(dat$attr$rCT.tx)))
+                                      (is.na(dat$attr$rCT.tx) | dat$attr$rCT.tx == 0)))
   idsUCT_prep_tx <- intersect(idsSTI_screen,
                               which(dat$attr$uCT == 1 &
                                       dat$attr$uCT.infTime < at &
-                                      is.na(dat$attr$uCT.tx)))
+                                      (is.na(dat$attr$uCT.tx) | dat$attr$uCT.tx == 0)))
 
   txRGC_prep <- idsRGC_prep_tx[which(rbinom(length(idsRGC_prep_tx), 1,
                                             prep.sti.prob.tx) == 1)]

From 41ed73dace13439f37f5ec203d86b48c4dba8318 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Thu, 8 Sep 2016 10:34:11 -0400
Subject: [PATCH 20/40] Update Rcpp exports again

---
 R/RcppExports.R     |  2 +-
 src/RcppExports.cpp | 18 +++++++++---------
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/R/RcppExports.R b/R/RcppExports.R
index 29b12d23..d3682776 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -1,4 +1,4 @@
-# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# This file was generated by Rcpp::compileAttributes
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 #' @title Aging Module
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index b5ca1d1a..7d2ad9b6 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -1,4 +1,4 @@
-// Generated by using Rcpp::compileAttributes() -> do not edit by hand
+// This file was generated by Rcpp::compileAttributes
 // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 #include <Rcpp.h>
@@ -9,23 +9,23 @@ using namespace Rcpp;
 List aging_msm(List dat, int at);
 RcppExport SEXP EpiModelHIV_aging_msm(SEXP datSEXP, SEXP atSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject rcpp_result_gen;
-    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::RObject __result;
+    Rcpp::RNGScope __rngScope;
     Rcpp::traits::input_parameter< List >::type dat(datSEXP);
     Rcpp::traits::input_parameter< int >::type at(atSEXP);
-    rcpp_result_gen = Rcpp::wrap(aging_msm(dat, at));
-    return rcpp_result_gen;
+    __result = Rcpp::wrap(aging_msm(dat, at));
+    return __result;
 END_RCPP
 }
 // aging_het
 List aging_het(List dat, int at);
 RcppExport SEXP EpiModelHIV_aging_het(SEXP datSEXP, SEXP atSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject rcpp_result_gen;
-    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::RObject __result;
+    Rcpp::RNGScope __rngScope;
     Rcpp::traits::input_parameter< List >::type dat(datSEXP);
     Rcpp::traits::input_parameter< int >::type at(atSEXP);
-    rcpp_result_gen = Rcpp::wrap(aging_het(dat, at));
-    return rcpp_result_gen;
+    __result = Rcpp::wrap(aging_het(dat, at));
+    return __result;
 END_RCPP
 }

From fecda9b66dd7ce9aca184db3ee55f20ef44fd1b4 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Thu, 8 Sep 2016 14:17:11 -0400
Subject: [PATCH 21/40] Final Rcpp doc updates

---
 R/RcppExports.R     |  2 +-
 src/RcppExports.cpp | 18 +++++++++---------
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/R/RcppExports.R b/R/RcppExports.R
index d3682776..29b12d23 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -1,4 +1,4 @@
-# This file was generated by Rcpp::compileAttributes
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 #' @title Aging Module
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 7d2ad9b6..b5ca1d1a 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -1,4 +1,4 @@
-// This file was generated by Rcpp::compileAttributes
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
 // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 #include <Rcpp.h>
@@ -9,23 +9,23 @@ using namespace Rcpp;
 List aging_msm(List dat, int at);
 RcppExport SEXP EpiModelHIV_aging_msm(SEXP datSEXP, SEXP atSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< List >::type dat(datSEXP);
     Rcpp::traits::input_parameter< int >::type at(atSEXP);
-    __result = Rcpp::wrap(aging_msm(dat, at));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(aging_msm(dat, at));
+    return rcpp_result_gen;
 END_RCPP
 }
 // aging_het
 List aging_het(List dat, int at);
 RcppExport SEXP EpiModelHIV_aging_het(SEXP datSEXP, SEXP atSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< List >::type dat(datSEXP);
     Rcpp::traits::input_parameter< int >::type at(atSEXP);
-    __result = Rcpp::wrap(aging_het(dat, at));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(aging_het(dat, at));
+    return rcpp_result_gen;
 END_RCPP
 }

From c8a4c0344f57a481975a99e149798d5c4d6fc4da Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Sun, 11 Sep 2016 17:01:13 -0400
Subject: [PATCH 22/40] Updating base parameters for calibrated model

---
 R/params.R       | 48 ++++++++++++++++++++++++------------------------
 man/init_msm.Rd  |  4 ++--
 man/param_msm.Rd | 29 ++++++++++++++---------------
 3 files changed, 40 insertions(+), 41 deletions(-)

diff --git a/R/params.R b/R/params.R
index 93c5ce79..8a735353 100644
--- a/R/params.R
+++ b/R/params.R
@@ -371,7 +371,7 @@ param_msm <- function(nwstats,
                       base.ai.pers.BB.rate = 0.11,
                       base.ai.pers.BW.rate = 0.16,
                       base.ai.pers.WW.rate = 0.14,
-                      ai.scale = 1,
+                      ai.scale = 1.15,
 
                       cond.main.BB.prob = 0.38,
                       cond.main.BW.prob = 0.10,
@@ -415,28 +415,28 @@ param_msm <- function(nwstats,
                       rcomp.main.only = FALSE,
                       rcomp.discl.only = FALSE,
 
-                      rgc.tprob = 0.45717564,
-                      ugc.tprob = 0.18828959,
-                      rct.tprob = 0.40582277,
-                      uct.tprob = 0.16310650,
+                      rgc.tprob = 0.3633513,
+                      ugc.tprob = 0.2479574,
+                      rct.tprob = 0.3123091,
+                      uct.tprob = 0.217685,
 
-                      rgc.sympt.prob = 0.07427077,
-                      ugc.sympt.prob = 0.86191181,
-                      rct.sympt.prob = 0.08108045,
-                      uct.sympt.prob = 0.91586663,
+                      rgc.sympt.prob = 0.06789649,
+                      ugc.sympt.prob = 0.8443316,
+                      rct.sympt.prob = 0.06979015,
+                      uct.sympt.prob = 0.9011038,
 
-                      rgc.dur.asympt = 42.82357437,
-                      ugc.dur.asympt = 33.43336473,
+                      rgc.dur.asympt = 37.76041,
+                      ugc.dur.asympt = 34.72606,
                       gc.dur.tx = 2,
                       gc.dur.ntx = NULL,
 
-                      rct.dur.asympt = 54.80423283,
-                      uct.dur.asympt = 46.20741625,
+                      rct.dur.asympt = 44.69662,
+                      uct.dur.asympt = 45.3775,
                       ct.dur.tx = 2,
                       ct.dur.ntx = NULL,
 
-                      gc.prob.cease = 0.03344859,
-                      ct.prob.cease = 0.03344859,
+                      gc.prob.cease = 0,
+                      ct.prob.cease = 0,
 
                       gc.sympt.prob.tx = 0.90,
                       ct.sympt.prob.tx = 0.85,
@@ -449,11 +449,11 @@ param_msm <- function(nwstats,
 
                       sti.cond.rr = 0.3,
 
-                      hiv.rgc.rr = 2.83399215,
-                      hiv.ugc.rr = 1.56832041,
-                      hiv.rct.rr = 2.83399215,
-                      hiv.uct.rr = 1.56832041,
-                      hiv.dual.rr = 0,
+                      hiv.rgc.rr = 2.878016,
+                      hiv.ugc.rr = 1.86111,
+                      hiv.rct.rr = 2.878016,
+                      hiv.uct.rr = 1.86111,
+                      hiv.dual.rr = 0.2,
                       ...) {
 
   p <- get_args(formal.args = formals(sys.function()),
@@ -571,10 +571,10 @@ param_msm <- function(nwstats,
 init_msm <- function(nwstats,
                      prev.B = 0.253,
                      prev.W = 0.253,
-                     prev.ugc = 0.01,
-                     prev.rgc = 0.05,
-                     prev.uct = 0.01,
-                     prev.rct = 0.05,
+                     prev.ugc = 0.005,
+                     prev.rgc = 0.005,
+                     prev.uct = 0.013,
+                     prev.rct = 0.013,
                      ...) {
 
   p <- get_args(formal.args = formals(sys.function()),
diff --git a/man/init_msm.Rd b/man/init_msm.Rd
index 8ba09a00..5b9eeef3 100644
--- a/man/init_msm.Rd
+++ b/man/init_msm.Rd
@@ -4,8 +4,8 @@
 \alias{init_msm}
 \title{Epidemic Model Initial Conditions}
 \usage{
-init_msm(nwstats, prev.B = 0.253, prev.W = 0.253, prev.ugc = 0.01,
-  prev.rgc = 0.05, prev.uct = 0.01, prev.rct = 0.05, ...)
+init_msm(nwstats, prev.B = 0.253, prev.W = 0.253, prev.ugc = 0.005,
+  prev.rgc = 0.005, prev.uct = 0.013, prev.rct = 0.013, ...)
 }
 \arguments{
 \item{nwstats}{Target statistics for the network model. An object of class
diff --git a/man/param_msm.Rd b/man/param_msm.Rd
index e6416968..f8aa6aa6 100644
--- a/man/param_msm.Rd
+++ b/man/param_msm.Rd
@@ -32,8 +32,8 @@ param_msm(nwstats, race.method = 1, last.neg.test.B.int = 301,
   ccr5.heteroz.rr = 0.3, num.inst.ai.classes = 1,
   base.ai.main.BB.rate = 0.17, base.ai.main.BW.rate = 0.26,
   base.ai.main.WW.rate = 0.23, base.ai.pers.BB.rate = 0.11,
-  base.ai.pers.BW.rate = 0.16, base.ai.pers.WW.rate = 0.14, ai.scale = 1,
-  cond.main.BB.prob = 0.38, cond.main.BW.prob = 0.1,
+  base.ai.pers.BW.rate = 0.16, base.ai.pers.WW.rate = 0.14,
+  ai.scale = 1.15, cond.main.BB.prob = 0.38, cond.main.BW.prob = 0.1,
   cond.main.WW.prob = 0.15, cond.pers.always.prob = 0.216,
   cond.pers.BB.prob = 0.26, cond.pers.BW.prob = 0.26,
   cond.pers.WW.prob = 0.26, cond.inst.always.prob = 0.326,
@@ -49,19 +49,18 @@ param_msm(nwstats, race.method = 1, last.neg.test.B.int = 301,
   prep.cov.method = "curr", prep.cov.rate = 1, prep.tst.int = 90,
   prep.risk.int = 182, prep.risk.reassess = TRUE, rcomp.prob = 0,
   rcomp.adh.groups = 0:3, rcomp.main.only = FALSE,
-  rcomp.discl.only = FALSE, rgc.tprob = 0.45717564,
-  ugc.tprob = 0.18828959, rct.tprob = 0.40582277, uct.tprob = 0.1631065,
-  rgc.sympt.prob = 0.07427077, ugc.sympt.prob = 0.86191181,
-  rct.sympt.prob = 0.08108045, uct.sympt.prob = 0.91586663,
-  rgc.dur.asympt = 42.82357437, ugc.dur.asympt = 33.43336473,
-  gc.dur.tx = 2, gc.dur.ntx = NULL, rct.dur.asympt = 54.80423283,
-  uct.dur.asympt = 46.20741625, ct.dur.tx = 2, ct.dur.ntx = NULL,
-  gc.prob.cease = 0.03344859, ct.prob.cease = 0.03344859,
-  gc.sympt.prob.tx = 0.9, ct.sympt.prob.tx = 0.85, gc.asympt.prob.tx = 0,
-  ct.asympt.prob.tx = 0, prep.sti.screen.int = 182, prep.sti.prob.tx = 1,
-  prep.continue.stand.tx = TRUE, sti.cond.rr = 0.3,
-  hiv.rgc.rr = 2.83399215, hiv.ugc.rr = 1.56832041,
-  hiv.rct.rr = 2.83399215, hiv.uct.rr = 1.56832041, hiv.dual.rr = 0, ...)
+  rcomp.discl.only = FALSE, rgc.tprob = 0.3633513, ugc.tprob = 0.2479574,
+  rct.tprob = 0.3123091, uct.tprob = 0.217685,
+  rgc.sympt.prob = 0.06789649, ugc.sympt.prob = 0.8443316,
+  rct.sympt.prob = 0.06979015, uct.sympt.prob = 0.9011038,
+  rgc.dur.asympt = 37.76041, ugc.dur.asympt = 34.72606, gc.dur.tx = 2,
+  gc.dur.ntx = NULL, rct.dur.asympt = 44.69662, uct.dur.asympt = 45.3775,
+  ct.dur.tx = 2, ct.dur.ntx = NULL, gc.prob.cease = 0,
+  ct.prob.cease = 0, gc.sympt.prob.tx = 0.9, ct.sympt.prob.tx = 0.85,
+  gc.asympt.prob.tx = 0, ct.asympt.prob.tx = 0, prep.sti.screen.int = 182,
+  prep.sti.prob.tx = 1, prep.continue.stand.tx = TRUE, sti.cond.rr = 0.3,
+  hiv.rgc.rr = 2.878016, hiv.ugc.rr = 1.86111, hiv.rct.rr = 2.878016,
+  hiv.uct.rr = 1.86111, hiv.dual.rr = 0.2, ...)
 }
 \arguments{
 \item{nwstats}{Target statistics for the network model. An object of class

From 989298edc8d66997712fbc7ee6919dc642297616 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Mon, 12 Sep 2016 09:57:43 -0400
Subject: [PATCH 23/40] Update .gitignore

---
 .gitignore | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitignore b/.gitignore
index 4dcad8d2..2a83a5bb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,4 @@ src/*.o
 src/*.so
 src/*.dll
 *.Rproj
+.DS_Store

From 0acd94775a4c3221f760dbb4a0e83b67a4b47842 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Wed, 14 Sep 2016 10:57:08 -0400
Subject: [PATCH 24/40] Add proportion of asympt and rect STIs treated

---
 R/mod.sti.R | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/R/mod.sti.R b/R/mod.sti.R
index e05d769d..61385b0a 100644
--- a/R/mod.sti.R
+++ b/R/mod.sti.R
@@ -596,6 +596,7 @@ sti_tx <- function(dat, at) {
   txUCT <- c(txUCT, txUCT_prep)
 
 
+
   # update attributes
   dat$attr$rGC.tx[c(idsRGC_tx, idsRGC_prep_tx)] <- 0
   dat$attr$rGC.tx[txRGC] <- 1
@@ -624,5 +625,29 @@ sti_tx <- function(dat, at) {
   dat$epi$txGC[at] <- length(txRGC) + length(txUGC)
   dat$epi$txCT[at] <- length(txRCT) + length(txUCT)
 
+  if (is.null(dat$epi$prop.GC.asympt.tx)) {
+    dat$epi$prop.GC.asympt.tx <- rep(NA, length(dat$epi$num))
+    dat$epi$prop.CT.asympt.tx <- rep(NA, length(dat$epi$num))
+    dat$epi$prop.rGC.tx <- rep(NA, length(dat$epi$num))
+    dat$epi$prop.rCT.tx <- rep(NA, length(dat$epi$num))
+  }
+  dat$epi$prop.GC.asympt.tx[at] <-
+    length(union(intersect(txRGC, which(dat$attr$rGC.sympt == 0)),
+                 intersect(txUGC, which(dat$attr$uGC.sympt == 0)))) /
+    length(union(union(idsRGC_tx_asympt,
+                       intersect(idsRGC_prep_tx, which(dat$attr$rGC.sympt == 0))),
+                 union(idsUGC_tx_asympt,
+                       intersect(idsUGC_prep_tx, which(dat$attr$uGC.sympt == 0)))))
+  dat$epi$prop.CT.asympt.tx[at] <-
+    length(union(intersect(txRCT, which(dat$attr$rCT.sympt == 0)),
+                 intersect(txUCT, which(dat$attr$uCT.sympt == 0)))) /
+    length(union(union(idsRCT_tx_asympt,
+                       intersect(idsRCT_prep_tx, which(dat$attr$rCT.sympt == 0))),
+                 union(idsUCT_tx_asympt,
+                       intersect(idsUCT_prep_tx, which(dat$attr$uCT.sympt == 0)))))
+
+  dat$epi$prop.rGC.tx[at] <- length(txRGC) / length(union(idsRGC_tx, idsRGC_prep_tx))
+  dat$epi$prop.rCT.tx[at] <- length(txRCT) / length(union(idsRCT_tx, idsRCT_prep_tx))
+
   return(dat)
 }

From 830e8b9beda79dfb53e3302e17d3911082c86296 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Wed, 14 Sep 2016 16:14:41 -0400
Subject: [PATCH 25/40] Fix bug with STI dx calcs

Change NaN to 0
---
 R/mod.sti.R | 18 ++++++++++++++----
 1 file changed, 14 insertions(+), 4 deletions(-)

diff --git a/R/mod.sti.R b/R/mod.sti.R
index 61385b0a..f518c398 100644
--- a/R/mod.sti.R
+++ b/R/mod.sti.R
@@ -631,23 +631,33 @@ sti_tx <- function(dat, at) {
     dat$epi$prop.rGC.tx <- rep(NA, length(dat$epi$num))
     dat$epi$prop.rCT.tx <- rep(NA, length(dat$epi$num))
   }
-  dat$epi$prop.GC.asympt.tx[at] <-
+  prop.GC.asympt.tx <-
     length(union(intersect(txRGC, which(dat$attr$rGC.sympt == 0)),
                  intersect(txUGC, which(dat$attr$uGC.sympt == 0)))) /
     length(union(union(idsRGC_tx_asympt,
                        intersect(idsRGC_prep_tx, which(dat$attr$rGC.sympt == 0))),
                  union(idsUGC_tx_asympt,
                        intersect(idsUGC_prep_tx, which(dat$attr$uGC.sympt == 0)))))
-  dat$epi$prop.CT.asympt.tx[at] <-
+  prop.GC.asympt.tx <- ifelse(is.nan(prop.GC.asympt.tx), 0, prop.GC.asympt.tx)
+  dat$epi$prop.GC.asympt.tx[at] <- prop.GC.asympt.tx
+
+  prop.CT.asympt.tx <-
     length(union(intersect(txRCT, which(dat$attr$rCT.sympt == 0)),
                  intersect(txUCT, which(dat$attr$uCT.sympt == 0)))) /
     length(union(union(idsRCT_tx_asympt,
                        intersect(idsRCT_prep_tx, which(dat$attr$rCT.sympt == 0))),
                  union(idsUCT_tx_asympt,
                        intersect(idsUCT_prep_tx, which(dat$attr$uCT.sympt == 0)))))
+  prop.CT.asympt.tx <- ifelse(is.nan(prop.CT.asympt.tx), 0, prop.CT.asympt.tx)
+  dat$epi$prop.CT.asympt.tx[at] <- prop.CT.asympt.tx
+
+  prop.rGC.tx <- length(txRGC) / length(union(idsRGC_tx, idsRGC_prep_tx))
+  prop.rGC.tx <- ifelse(is.nan(prop.rGC.tx), 0, prop.rGC.tx)
+  dat$epi$prop.rGC.tx[at] <- prop.rGC.tx
 
-  dat$epi$prop.rGC.tx[at] <- length(txRGC) / length(union(idsRGC_tx, idsRGC_prep_tx))
-  dat$epi$prop.rCT.tx[at] <- length(txRCT) / length(union(idsRCT_tx, idsRCT_prep_tx))
+  prop.rCT.tx <- length(txRCT) / length(union(idsRCT_tx, idsRCT_prep_tx))
+  prop.rCT.tx <- ifelse(is.nan(prop.rCT.tx), 0, prop.rCT.tx)
+  dat$epi$prop.rCT.tx[at] <- prop.rCT.tx
 
   return(dat)
 }

From 8bd804405342d0fb3c8d1de7ac69113082e7129b Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Tue, 27 Sep 2016 09:15:48 -0400
Subject: [PATCH 26/40] Update STI parameter values

---
 R/params.R       | 36 ++++++++++++++++++------------------
 man/param_msm.Rd | 24 ++++++++++++------------
 2 files changed, 30 insertions(+), 30 deletions(-)

diff --git a/R/params.R b/R/params.R
index 8a735353..38bcbe04 100644
--- a/R/params.R
+++ b/R/params.R
@@ -415,23 +415,23 @@ param_msm <- function(nwstats,
                       rcomp.main.only = FALSE,
                       rcomp.discl.only = FALSE,
 
-                      rgc.tprob = 0.3633513,
-                      ugc.tprob = 0.2479574,
-                      rct.tprob = 0.3123091,
-                      uct.tprob = 0.217685,
-
-                      rgc.sympt.prob = 0.06789649,
-                      ugc.sympt.prob = 0.8443316,
-                      rct.sympt.prob = 0.06979015,
-                      uct.sympt.prob = 0.9011038,
-
-                      rgc.dur.asympt = 37.76041,
-                      ugc.dur.asympt = 34.72606,
+                      rgc.tprob = 0.357698,
+                      ugc.tprob = 0.248095,
+                      rct.tprob = 0.321597,
+                      uct.tprob = 0.212965,
+
+                      rgc.sympt.prob = 0.076975,
+                      ugc.sympt.prob = 0.824368,
+                      rct.sympt.prob = 0.103517,
+                      uct.sympt.prob = 0.885045,
+
+                      rgc.dur.asympt = 34.93723,
+                      ugc.dur.asympt = 36.48328,
                       gc.dur.tx = 2,
                       gc.dur.ntx = NULL,
 
-                      rct.dur.asympt = 44.69662,
-                      uct.dur.asympt = 45.3775,
+                      rct.dur.asympt = 45.02343,
+                      uct.dur.asympt = 44.88562,
                       ct.dur.tx = 2,
                       ct.dur.ntx = NULL,
 
@@ -449,10 +449,10 @@ param_msm <- function(nwstats,
 
                       sti.cond.rr = 0.3,
 
-                      hiv.rgc.rr = 2.878016,
-                      hiv.ugc.rr = 1.86111,
-                      hiv.rct.rr = 2.878016,
-                      hiv.uct.rr = 1.86111,
+                      hiv.rgc.rr = 2.644584,
+                      hiv.ugc.rr = 1.69434,
+                      hiv.rct.rr = 2.644584,
+                      hiv.uct.rr = 1.69434,
                       hiv.dual.rr = 0.2,
                       ...) {
 
diff --git a/man/param_msm.Rd b/man/param_msm.Rd
index f8aa6aa6..2c8d8050 100644
--- a/man/param_msm.Rd
+++ b/man/param_msm.Rd
@@ -49,18 +49,18 @@ param_msm(nwstats, race.method = 1, last.neg.test.B.int = 301,
   prep.cov.method = "curr", prep.cov.rate = 1, prep.tst.int = 90,
   prep.risk.int = 182, prep.risk.reassess = TRUE, rcomp.prob = 0,
   rcomp.adh.groups = 0:3, rcomp.main.only = FALSE,
-  rcomp.discl.only = FALSE, rgc.tprob = 0.3633513, ugc.tprob = 0.2479574,
-  rct.tprob = 0.3123091, uct.tprob = 0.217685,
-  rgc.sympt.prob = 0.06789649, ugc.sympt.prob = 0.8443316,
-  rct.sympt.prob = 0.06979015, uct.sympt.prob = 0.9011038,
-  rgc.dur.asympt = 37.76041, ugc.dur.asympt = 34.72606, gc.dur.tx = 2,
-  gc.dur.ntx = NULL, rct.dur.asympt = 44.69662, uct.dur.asympt = 45.3775,
-  ct.dur.tx = 2, ct.dur.ntx = NULL, gc.prob.cease = 0,
-  ct.prob.cease = 0, gc.sympt.prob.tx = 0.9, ct.sympt.prob.tx = 0.85,
-  gc.asympt.prob.tx = 0, ct.asympt.prob.tx = 0, prep.sti.screen.int = 182,
-  prep.sti.prob.tx = 1, prep.continue.stand.tx = TRUE, sti.cond.rr = 0.3,
-  hiv.rgc.rr = 2.878016, hiv.ugc.rr = 1.86111, hiv.rct.rr = 2.878016,
-  hiv.uct.rr = 1.86111, hiv.dual.rr = 0.2, ...)
+  rcomp.discl.only = FALSE, rgc.tprob = 0.357698, ugc.tprob = 0.248095,
+  rct.tprob = 0.321597, uct.tprob = 0.212965, rgc.sympt.prob = 0.076975,
+  ugc.sympt.prob = 0.824368, rct.sympt.prob = 0.103517,
+  uct.sympt.prob = 0.885045, rgc.dur.asympt = 34.93723,
+  ugc.dur.asympt = 36.48328, gc.dur.tx = 2, gc.dur.ntx = NULL,
+  rct.dur.asympt = 45.02343, uct.dur.asympt = 44.88562, ct.dur.tx = 2,
+  ct.dur.ntx = NULL, gc.prob.cease = 0, ct.prob.cease = 0,
+  gc.sympt.prob.tx = 0.9, ct.sympt.prob.tx = 0.85, gc.asympt.prob.tx = 0,
+  ct.asympt.prob.tx = 0, prep.sti.screen.int = 182, prep.sti.prob.tx = 1,
+  prep.continue.stand.tx = TRUE, sti.cond.rr = 0.3, hiv.rgc.rr = 2.644584,
+  hiv.ugc.rr = 1.69434, hiv.rct.rr = 2.644584, hiv.uct.rr = 1.69434,
+  hiv.dual.rr = 0.2, ...)
 }
 \arguments{
 \item{nwstats}{Target statistics for the network model. An object of class

From b91b4c82dde6bf7eb5e548ee655d6855ca4da3c2 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Tue, 27 Sep 2016 13:25:13 -0400
Subject: [PATCH 27/40] Add prep specific tx attr

---
 R/mod.sti.R | 166 ++++++++++++++++++++++++++++++++--------------------
 1 file changed, 101 insertions(+), 65 deletions(-)

diff --git a/R/mod.sti.R b/R/mod.sti.R
index f518c398..c83d88bd 100644
--- a/R/mod.sti.R
+++ b/R/mod.sti.R
@@ -296,7 +296,8 @@ sti_trans <- function(dat, at) {
 #'
 sti_recov <- function(dat, at) {
 
-  # Parameters
+  # Parameters ----------------------------------------------------------
+
   rgc.dur.asympt <- dat$param$rgc.dur.asympt
   ugc.dur.asympt <- dat$param$ugc.dur.asympt
   gc.dur.tx <- dat$param$gc.dur.tx
@@ -308,47 +309,58 @@ sti_recov <- function(dat, at) {
   ct.dur.ntx <- dat$param$ct.dur.ntx
 
 
-  # GC recovery
+  # GC Recovery ---------------------------------------------------------
+
+  # Asymptomatic untreated
   idsRGC_asympt_ntx <- which(dat$attr$rGC == 1 &
                              dat$attr$rGC.infTime < at &
                              dat$attr$rGC.sympt == 0 &
-                             (is.na(dat$attr$rGC.tx) | dat$attr$rGC.tx == 0))
+                             (is.na(dat$attr$rGC.tx) | dat$attr$rGC.tx == 0) &
+                             (is.na(dat$attr$rGC.tx.prep) | dat$attr$rGC.tx.prep == 0))
   idsUGC_asympt_ntx <- which(dat$attr$uGC == 1 &
                              dat$attr$uGC.infTime < at &
                              dat$attr$uGC.sympt == 0 &
-                             (is.na(dat$attr$uGC.tx) | dat$attr$uGC.tx == 0))
+                             (is.na(dat$attr$uGC.tx) | dat$attr$uGC.tx == 0) &
+                             (is.na(dat$attr$uGC.tx.prep) | dat$attr$uGC.tx.prep == 0))
+
+  recovRGC_asympt_ntx <- idsRGC_asympt_ntx[which(rbinom(length(idsRGC_asympt_ntx), 1,
+                                                        1/rgc.dur.asympt) == 1)]
+  recovUGC_asympt_ntx <- idsUGC_asympt_ntx[which(rbinom(length(idsUGC_asympt_ntx), 1,
+                                                        1/ugc.dur.asympt) == 1)]
+
+  # Symptomatic untreated
   idsRGC_sympt_ntx <- which(dat$attr$rGC == 1 &
                             dat$attr$rGC.infTime < at &
                             dat$attr$rGC.sympt == 1 &
-                            (is.na(dat$attr$rGC.tx) | dat$attr$rGC.tx == 0))
+                            (is.na(dat$attr$rGC.tx) | dat$attr$rGC.tx == 0) &
+                            (is.na(dat$attr$rGC.tx.prep) | dat$attr$rGC.tx.prep == 0))
   idsUGC_sympt_ntx <- which(dat$attr$uGC == 1 &
                             dat$attr$uGC.infTime < at &
                             dat$attr$uGC.sympt == 1 &
-                            (is.na(dat$attr$uGC.tx) | dat$attr$uGC.tx == 0))
-  idsRGC_tx <- which(dat$attr$rGC == 1 &
-                     dat$attr$rGC.infTime < at &
-                     dat$attr$rGC.tx == 1)
-  idsUGC_tx <- which(dat$attr$uGC == 1 &
-                     dat$attr$uGC.infTime < at &
-                     dat$attr$uGC.tx == 1)
-
-  recovRGC_asympt_ntx <- idsRGC_asympt_ntx[which(rbinom(length(idsRGC_asympt_ntx), 1,
-                                                 1/rgc.dur.asympt) == 1)]
-  recovUGC_asympt_ntx <- idsUGC_asympt_ntx[which(rbinom(length(idsUGC_asympt_ntx), 1,
-                                                 1/ugc.dur.asympt) == 1)]
+                            (is.na(dat$attr$uGC.tx) | dat$attr$uGC.tx == 0) &
+                            (is.na(dat$attr$uGC.tx.prep) | dat$attr$uGC.tx.prep == 0))
 
+  # If parameter is null, uses recovery rate of asytomatic untreated
   if (!is.null(gc.dur.ntx)) {
     recovRGC_sympt_ntx <- idsRGC_sympt_ntx[which(rbinom(length(idsRGC_sympt_ntx), 1,
-                                                 1/gc.dur.ntx) == 1)]
+                                                        1/gc.dur.ntx) == 1)]
     recovUGC_sympt_ntx <- idsUGC_sympt_ntx[which(rbinom(length(idsUGC_sympt_ntx), 1,
-                                                 1/gc.dur.ntx) == 1)]
+                                                        1/gc.dur.ntx) == 1)]
   } else {
     recovRGC_sympt_ntx <- idsRGC_sympt_ntx[which(rbinom(length(idsRGC_sympt_ntx), 1,
-                                                 1/rgc.dur.asympt) == 1)]
+                                                        1/rgc.dur.asympt) == 1)]
     recovUGC_sympt_ntx <- idsUGC_sympt_ntx[which(rbinom(length(idsUGC_sympt_ntx), 1,
-                                                 1/ugc.dur.asympt) == 1)]
+                                                        1/ugc.dur.asympt) == 1)]
   }
 
+  # Treated (asymptomatic and symptomatic)
+  idsRGC_tx <- which(dat$attr$rGC == 1 &
+                     dat$attr$rGC.infTime < at &
+                     (dat$attr$rGC.tx == 1 | dat$attr$rGC.tx.prep == 1))
+  idsUGC_tx <- which(dat$attr$uGC == 1 &
+                     dat$attr$uGC.infTime < at &
+                     (dat$attr$uGC.tx == 1 | dat$attr$uGC.tx.prep == 1))
+
   recovRGC_tx <- idsRGC_tx[which(rbinom(length(idsRGC_tx), 1,
                                         1/gc.dur.tx) == 1)]
   recovUGC_tx <- idsUGC_tx[which(rbinom(length(idsUGC_tx), 1,
@@ -361,60 +373,75 @@ sti_recov <- function(dat, at) {
   dat$attr$rGC.sympt[recovRGC] <- NA
   dat$attr$rGC.infTime[recovRGC] <- NA
   dat$attr$rGC.tx[recovRGC] <- NA
+  dat$attr$rGC.tx.prep[recovRGC] <- NA
 
   dat$attr$uGC[recovUGC] <- 0
   dat$attr$uGC.sympt[recovUGC] <- NA
   dat$attr$uGC.infTime[recovUGC] <- NA
   dat$attr$uGC.tx[recovUGC] <- NA
+  dat$attr$uGC.tx.prep[recovUGC] <- NA
 
   dat$attr$GC.cease[c(recovRGC, recovUGC)] <- NA
 
-  # CT recovery
+
+
+  # CT Recovery ---------------------------------------------------------
+
+  # Asymptomatic untreated
   idsRCT_asympt_ntx <- which(dat$attr$rCT == 1 &
                              dat$attr$rCT.infTime < at &
                              dat$attr$rCT.sympt == 0 &
-                             (is.na(dat$attr$rCT.tx) | dat$attr$rCT.tx == 0))
+                             (is.na(dat$attr$rCT.tx) | dat$attr$rCT.tx == 0) &
+                             (is.na(dat$attr$rCT.tx.prep) | dat$attr$rCT.tx.prep == 0))
   idsUCT_asympt_ntx <- which(dat$attr$uCT == 1 &
                              dat$attr$uCT.infTime < at &
                              dat$attr$uCT.sympt == 0 &
-                             (is.na(dat$attr$uCT.tx) | dat$attr$uCT.tx == 0))
+                             (is.na(dat$attr$uCT.tx) | dat$attr$uCT.tx == 0) &
+                             (is.na(dat$attr$uCT.tx.prep) | dat$attr$uCT.tx.prep == 0))
+
+  recovRCT_asympt_ntx <- idsRCT_asympt_ntx[which(rbinom(length(idsRCT_asympt_ntx),
+                                                        1, 1/rct.dur.asympt) == 1)]
+  recovUCT_asympt_ntx <- idsUCT_asympt_ntx[which(rbinom(length(idsUCT_asympt_ntx),
+                                                        1, 1/uct.dur.asympt) == 1)]
+
+  # Symptomatic untreated
   idsRCT_sympt_ntx <- which(dat$attr$rCT == 1 &
                             dat$attr$rCT.infTime < at &
                             dat$attr$rCT.sympt == 1 &
-                            (is.na(dat$attr$rCT.tx) | dat$attr$rCT.tx == 0))
+                            (is.na(dat$attr$rCT.tx) | dat$attr$rCT.tx == 0) &
+                            (is.na(dat$attr$rCT.tx.prep) | dat$attr$rCT.tx.prep == 0))
   idsUCT_sympt_ntx <- which(dat$attr$uCT == 1 &
                             dat$attr$uCT.infTime < at &
                             dat$attr$uCT.sympt == 1 &
-                            (is.na(dat$attr$uCT.tx) | dat$attr$uCT.tx == 0))
-  idsRCT_tx <- which(dat$attr$rCT == 1 &
-                     dat$attr$rCT.infTime < at &
-                     dat$attr$rCT.tx == 1)
-  idsUCT_tx <- which(dat$attr$uCT == 1 &
-                     dat$attr$uCT.infTime < at &
-                     dat$attr$uCT.tx == 1)
-
-  recovRCT_asympt_ntx <- idsRCT_asympt_ntx[which(rbinom(length(idsRCT_asympt_ntx),
-                                                 1, 1/rct.dur.asympt) == 1)]
-  recovUCT_asympt_ntx <- idsUCT_asympt_ntx[which(rbinom(length(idsUCT_asympt_ntx),
-                                                 1, 1/uct.dur.asympt) == 1)]
+                            (is.na(dat$attr$uCT.tx) | dat$attr$uCT.tx == 0) &
+                            (is.na(dat$attr$uCT.tx.prep) | dat$attr$uCT.tx.prep == 0))
 
   if (!is.null(ct.dur.ntx)) {
     recovRCT_sympt_ntx <- idsRCT_sympt_ntx[which(rbinom(length(idsRCT_sympt_ntx),
-                                            1, 1/ct.dur.ntx) == 1)]
+                                                        1, 1/ct.dur.ntx) == 1)]
     recovUCT_sympt_ntx <- idsUCT_sympt_ntx[which(rbinom(length(idsUCT_sympt_ntx),
-                                            1, 1/ct.dur.ntx) == 1)]
+                                                        1, 1/ct.dur.ntx) == 1)]
   } else {
     recovRCT_sympt_ntx <- idsRCT_sympt_ntx[which(rbinom(length(idsRCT_sympt_ntx),
-                                            1, 1/rct.dur.asympt) == 1)]
+                                                        1, 1/rct.dur.asympt) == 1)]
     recovUCT_sympt_ntx <- idsUCT_sympt_ntx[which(rbinom(length(idsUCT_sympt_ntx),
-                                            1, 1/uct.dur.asympt) == 1)]
+                                                        1, 1/uct.dur.asympt) == 1)]
   }
 
+  # Treated (asymptomatic and symptomatic)
+  idsRCT_tx <- which(dat$attr$rCT == 1 &
+                     dat$attr$rCT.infTime < at &
+                     (dat$attr$rCT.tx == 1 | dat$attr$rCT.tx.prep == 1))
+  idsUCT_tx <- which(dat$attr$uCT == 1 &
+                     dat$attr$uCT.infTime < at &
+                     (dat$attr$uCT.tx == 1 | dat$attr$uCT.tx.prep == 1))
+
   recovRCT_tx <- idsRCT_tx[which(rbinom(length(idsRCT_tx),
                                         1, 1/ct.dur.tx) == 1)]
   recovUCT_tx <- idsUCT_tx[which(rbinom(length(idsUCT_tx),
                                         1, 1/ct.dur.tx) == 1)]
 
+
   recovRCT <- c(recovRCT_asympt_ntx, recovRCT_sympt_ntx, recovRCT_tx)
   recovUCT <- c(recovUCT_asympt_ntx, recovUCT_sympt_ntx, recovUCT_tx)
 
@@ -422,19 +449,21 @@ sti_recov <- function(dat, at) {
   dat$attr$rCT.sympt[recovRCT] <- NA
   dat$attr$rCT.infTime[recovRCT] <- NA
   dat$attr$rCT.tx[recovRCT] <- NA
+  dat$attr$rCT.tx.prep[recovRCT] <- NA
 
   dat$attr$uCT[recovUCT] <- 0
   dat$attr$uCT.sympt[recovUCT] <- NA
   dat$attr$uCT.infTime[recovUCT] <- NA
   dat$attr$uCT.tx[recovUCT] <- NA
+  dat$attr$uCT.tx.prep[recovUCT] <- NA
 
   dat$attr$CT.cease[c(recovRCT, recovUCT)] <- NA
 
   # Summary stats
-  dat$epi$recov.rgc[at] <- length(recovRGC)
-  dat$epi$recov.ugc[at] <- length(recovUGC)
-  dat$epi$recov.rct[at] <- length(recovRCT)
-  dat$epi$recov.uct[at] <- length(recovUCT)
+  dat$epi$recov.rgc[at] <- length(unique(recovRGC))
+  dat$epi$recov.ugc[at] <- length(unique(recovUGC))
+  dat$epi$recov.rct[at] <- length(unique(recovRCT))
+  dat$epi$recov.uct[at] <- length(unique(recovUCT))
 
   return(dat)
 }
@@ -562,24 +591,23 @@ sti_tx <- function(dat, at) {
 
   dat$attr$prepLastStiScreen[idsSTI_screen] <- at
 
-  ## TODO: current approach does not allow for fixed prep-specific non-treatment
-  ##       attribute, where the prep.sti.prob.tx parameter < 1 (default)
+
   idsRGC_prep_tx <- intersect(idsSTI_screen,
                               which(dat$attr$rGC == 1 &
                                       dat$attr$rGC.infTime < at &
-                                      (is.na(dat$attr$rGC.tx) | dat$attr$rGC.tx == 0)))
+                                      is.na(dat$attr$rGC.tx.prep)))
   idsUGC_prep_tx <- intersect(idsSTI_screen,
                               which(dat$attr$uGC == 1 &
                                       dat$attr$uGC.infTime < at &
-                                      (is.na(dat$attr$uGC.tx) | dat$attr$uGC.tx == 0)))
+                                      is.na(dat$attr$uGC.tx.prep)))
   idsRCT_prep_tx <- intersect(idsSTI_screen,
                               which(dat$attr$rCT == 1 &
                                       dat$attr$rCT.infTime < at &
-                                      (is.na(dat$attr$rCT.tx) | dat$attr$rCT.tx == 0)))
+                                      is.na(dat$attr$rCT.tx.prep)))
   idsUCT_prep_tx <- intersect(idsSTI_screen,
                               which(dat$attr$uCT == 1 &
                                       dat$attr$uCT.infTime < at &
-                                      (is.na(dat$attr$uCT.tx) | dat$attr$uCT.tx == 0)))
+                                      is.na(dat$attr$uCT.tx.prep)))
 
   txRGC_prep <- idsRGC_prep_tx[which(rbinom(length(idsRGC_prep_tx), 1,
                                             prep.sti.prob.tx) == 1)]
@@ -590,32 +618,40 @@ sti_tx <- function(dat, at) {
   txUCT_prep <- idsUCT_prep_tx[which(rbinom(length(idsUCT_prep_tx), 1,
                                             prep.sti.prob.tx) == 1)]
 
-  txRGC <- c(txRGC, txRGC_prep)
-  txUGC <- c(txUGC, txUGC_prep)
-  txRCT <- c(txRCT, txRCT_prep)
-  txUCT <- c(txUCT, txUCT_prep)
-
-
 
   # update attributes
-  dat$attr$rGC.tx[c(idsRGC_tx, idsRGC_prep_tx)] <- 0
+  dat$attr$rGC.tx[idsRGC_tx] <- 0
   dat$attr$rGC.tx[txRGC] <- 1
 
-  dat$attr$uGC.tx[c(idsUGC_tx, idsUGC_prep_tx)] <- 0
+  dat$attr$uGC.tx[idsUGC_tx] <- 0
   dat$attr$uGC.tx[txUGC] <- 1
 
-  dat$attr$rCT.tx[c(idsRCT_tx, idsRCT_prep_tx)] <- 0
+  dat$attr$rCT.tx[idsRCT_tx] <- 0
   dat$attr$rCT.tx[txRCT] <- 1
 
-  dat$attr$uCT.tx[c(idsUCT_tx, idsUCT_prep_tx)] <- 0
+  dat$attr$uCT.tx[idsUCT_tx] <- 0
   dat$attr$uCT.tx[txUCT] <- 1
 
+  dat$attr$rGC.tx.prep[idsRGC_prep_tx] <- 0
+  dat$attr$rGC.tx.prep[txRGC_prep] <- 1
+
+  dat$attr$uGC.tx.prep[idsUGC_prep_tx] <- 0
+  dat$attr$uGC.tx.prep[txUGC_prep] <- 1
+
+  dat$attr$rCT.tx.prep[idsRCT_prep_tx] <- 0
+  dat$attr$rCT.tx.prep[txRCT_prep] <- 1
+
+  dat$attr$uCT.tx.prep[idsUCT_prep_tx] <- 0
+  dat$attr$uCT.tx.prep[txUCT_prep] <- 1
+
+
   # add tx at other site
-  dat$attr$rGC.tx[which(dat$attr$uGC.tx == 1 & dat$attr$rGC == 1)] <- 1
-  dat$attr$uGC.tx[which(dat$attr$rGC.tx == 1 & dat$attr$uGC == 1)] <- 1
+  dat$attr$rGC.tx[which((dat$attr$uGC.tx == 1 | dat$attr$uGC.tx.prep == 1) & dat$attr$rGC == 1)] <- 1
+  dat$attr$uGC.tx[which((dat$attr$rGC.tx == 1 | dat$attr$rGC.tx.prep == 1) & dat$attr$uGC == 1)] <- 1
+
+  dat$attr$rCT.tx[which((dat$attr$uCT.tx == 1 | dat$attr$uCT.tx.prep == 1) & dat$attr$rCT == 1)] <- 1
+  dat$attr$uCT.tx[which((dat$attr$rCT.tx == 1 | dat$attr$rCT.tx.prep == 1) & dat$attr$uCT == 1)] <- 1
 
-  dat$attr$rCT.tx[which(dat$attr$uCT.tx == 1 & dat$attr$rCT == 1)] <- 1
-  dat$attr$uCT.tx[which(dat$attr$rCT.tx == 1 & dat$attr$uCT == 1)] <- 1
 
   # summary stats
   if (is.null(dat$epi$txGC)) {

From 2fb26e4bcf11ad0830761a1f32569eb3bc1128af Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Tue, 27 Sep 2016 13:25:32 -0400
Subject: [PATCH 28/40] New tx-related summary stats for paper

---
 R/mod.initialize.R |  3 ++
 R/mod.prevalence.R |  8 ++++
 R/mod.sti.R        | 91 ++++++++++++++++++++++++++--------------------
 3 files changed, 63 insertions(+), 39 deletions(-)

diff --git a/R/mod.initialize.R b/R/mod.initialize.R
index 7010216b..17309698 100644
--- a/R/mod.initialize.R
+++ b/R/mod.initialize.R
@@ -170,6 +170,7 @@ initialize_msm <- function(x, param, init, control, s) {
   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
@@ -196,8 +197,10 @@ initialize_msm <- function(x, param, init, control, s) {
   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)
 
diff --git a/R/mod.prevalence.R b/R/mod.prevalence.R
index 51d02292..350de2a6 100644
--- a/R/mod.prevalence.R
+++ b/R/mod.prevalence.R
@@ -96,6 +96,9 @@ prevalence_msm <- function(dat, at) {
     dat$epi$ir100.uct <- rNA
     dat$epi$ir100.ct <- rNA
 
+    dat$epi$ir100.sti <- rNA
+    dat$epi$incid.gcct.prep <- rNA
+
     dat$epi$recov.rgc <- rNA
     dat$epi$recov.ugc <- rNA
     dat$epi$recov.rct <- rNA
@@ -168,6 +171,11 @@ prevalence_msm <- function(dat, at) {
                                  sum(rCT == 0, na.rm = TRUE) +
                                  sum(uCT == 0, na.rm = TRUE))) * 5200
 
+  dat$epi$ir100.sti.prep[at] <- (dat$epi$incid.gcct.prep[at] /
+                                  (sum(rGC == 0 & prepStat == 1, na.rm = TRUE) +
+                                   sum(uGC == 0 & prepStat == 1, na.rm = TRUE) +
+                                   sum(rCT == 0 & prepStat == 1, na.rm = TRUE) +
+                                   sum(uCT == 0 & prepStat == 1, na.rm = TRUE))) * 5200
 
   return(dat)
 }
diff --git a/R/mod.sti.R b/R/mod.sti.R
index c83d88bd..f8154b54 100644
--- a/R/mod.sti.R
+++ b/R/mod.sti.R
@@ -259,6 +259,10 @@ sti_trans <- function(dat, at) {
   dat$epi$incid.uct[at] <- length(idsInf_uct)
   dat$epi$incid.ct[at] <- length(idsInf_rct) + length(idsInf_uct)
 
+  dat$epi$incid.gcct.prep[at] <- length(intersect(unique(c(idsInf_rgc, idsInf_ugc,
+                                                    idsInf_rct, idsInf_uct)),
+                                           which(dat$attr$prepStat == 1)))
+
   # Check all infected have all STI attributes
   stopifnot(all(!is.na(rGC.infTime[rGC == 1])),
             all(!is.na(rGC.sympt[rGC == 1])),
@@ -652,48 +656,57 @@ sti_tx <- function(dat, at) {
   dat$attr$rCT.tx[which((dat$attr$uCT.tx == 1 | dat$attr$uCT.tx.prep == 1) & dat$attr$rCT == 1)] <- 1
   dat$attr$uCT.tx[which((dat$attr$rCT.tx == 1 | dat$attr$rCT.tx.prep == 1) & dat$attr$uCT == 1)] <- 1
 
+  txRGC_all <- union(txRGC, txRGC_prep)
+  txUGC_all <- union(txUGC, txUGC_prep)
+  txRCT_all <- union(txRCT, txRCT_prep)
+  txUCT_all <- union(txUCT, txUCT_prep)
+
 
   # summary stats
-  if (is.null(dat$epi$txGC)) {
-    dat$epi$txGC <- rep(NA, length(dat$epi$num))
-    dat$epi$txCT <- rep(NA, length(dat$epi$num))
+  if (is.null(dat$epi$num.asympt.tx)) {
+    dat$epi$num.asympt.tx <- rep(NA, length(dat$epi$num))
+    dat$epi$num.asympt.cases <- rep(NA, length(dat$epi$num))
+    dat$epi$num.asympt.tx.prep <- rep(NA, length(dat$epi$num))
+    dat$epi$num.asympt.cases.prep <- rep(NA, length(dat$epi$num))
+    dat$epi$num.rect.tx <- rep(NA, length(dat$epi$num))
+    dat$epi$num.rect.cases <- rep(NA, length(dat$epi$num))
+    dat$epi$num.rect.tx.prep <- rep(NA, length(dat$epi$num))
+    dat$epi$num.rect.cases.prep <- rep(NA, length(dat$epi$num))
   }
-  dat$epi$txGC[at] <- length(txRGC) + length(txUGC)
-  dat$epi$txCT[at] <- length(txRCT) + length(txUCT)
-
-  if (is.null(dat$epi$prop.GC.asympt.tx)) {
-    dat$epi$prop.GC.asympt.tx <- rep(NA, length(dat$epi$num))
-    dat$epi$prop.CT.asympt.tx <- rep(NA, length(dat$epi$num))
-    dat$epi$prop.rGC.tx <- rep(NA, length(dat$epi$num))
-    dat$epi$prop.rCT.tx <- rep(NA, length(dat$epi$num))
-  }
-  prop.GC.asympt.tx <-
-    length(union(intersect(txRGC, which(dat$attr$rGC.sympt == 0)),
-                 intersect(txUGC, which(dat$attr$uGC.sympt == 0)))) /
-    length(union(union(idsRGC_tx_asympt,
-                       intersect(idsRGC_prep_tx, which(dat$attr$rGC.sympt == 0))),
-                 union(idsUGC_tx_asympt,
-                       intersect(idsUGC_prep_tx, which(dat$attr$uGC.sympt == 0)))))
-  prop.GC.asympt.tx <- ifelse(is.nan(prop.GC.asympt.tx), 0, prop.GC.asympt.tx)
-  dat$epi$prop.GC.asympt.tx[at] <- prop.GC.asympt.tx
-
-  prop.CT.asympt.tx <-
-    length(union(intersect(txRCT, which(dat$attr$rCT.sympt == 0)),
-                 intersect(txUCT, which(dat$attr$uCT.sympt == 0)))) /
-    length(union(union(idsRCT_tx_asympt,
-                       intersect(idsRCT_prep_tx, which(dat$attr$rCT.sympt == 0))),
-                 union(idsUCT_tx_asympt,
-                       intersect(idsUCT_prep_tx, which(dat$attr$uCT.sympt == 0)))))
-  prop.CT.asympt.tx <- ifelse(is.nan(prop.CT.asympt.tx), 0, prop.CT.asympt.tx)
-  dat$epi$prop.CT.asympt.tx[at] <- prop.CT.asympt.tx
-
-  prop.rGC.tx <- length(txRGC) / length(union(idsRGC_tx, idsRGC_prep_tx))
-  prop.rGC.tx <- ifelse(is.nan(prop.rGC.tx), 0, prop.rGC.tx)
-  dat$epi$prop.rGC.tx[at] <- prop.rGC.tx
-
-  prop.rCT.tx <- length(txRCT) / length(union(idsRCT_tx, idsRCT_prep_tx))
-  prop.rCT.tx <- ifelse(is.nan(prop.rCT.tx), 0, prop.rCT.tx)
-  dat$epi$prop.rCT.tx[at] <- prop.rCT.tx
+
+  asympt.tx <- c(intersect(txRGC_all, which(dat$attr$rGC.sympt == 0)),
+                 intersect(txUGC_all, which(dat$attr$uGC.sympt == 0)),
+                 intersect(txRCT_all, which(dat$attr$rCT.sympt == 0)),
+                 intersect(txUCT_all, which(dat$attr$uCT.sympt == 0)))
+  dat$epi$num.asympt.tx[at] <- length(unique(asympt.tx))
+  asympt.cases <- c(idsRGC_tx_asympt, intersect(idsRGC_prep_tx, which(dat$attr$rGC.sympt == 0)),
+                    idsUGC_tx_asympt, intersect(idsUGC_prep_tx, which(dat$attr$uGC.sympt == 0)),
+                    idsRCT_tx_asympt, intersect(idsRCT_prep_tx, which(dat$attr$rCT.sympt == 0)),
+                    idsUCT_tx_asympt, intersect(idsUCT_prep_tx, which(dat$attr$uCT.sympt == 0)))
+  dat$epi$num.asympt.cases[at] <- length(unique(asympt.cases))
+
+
+  asympt.tx.prep <- c(intersect(txRGC_prep, which(dat$attr$rGC.sympt == 0)),
+                      intersect(txUGC_prep, which(dat$attr$uGC.sympt == 0)),
+                      intersect(txRCT_prep, which(dat$attr$rCT.sympt == 0)),
+                      intersect(txUCT_prep, which(dat$attr$uCT.sympt == 0)))
+  dat$epi$num.asympt.tx.prep[at] <- length(unique(asympt.tx.prep))
+  asympt.cases.prep <- c(intersect(idsRGC_prep_tx, which(dat$attr$rGC.sympt == 0)),
+                         intersect(idsUGC_prep_tx, which(dat$attr$uGC.sympt == 0)),
+                         intersect(idsRCT_prep_tx, which(dat$attr$rCT.sympt == 0)),
+                         intersect(idsUCT_prep_tx, which(dat$attr$uCT.sympt == 0)))
+  dat$epi$num.asympt.cases.prep[at] <- length(unique(asympt.cases.prep))
+
+
+  rect.tx <- c(txRGC_all, txRCT_all)
+  dat$epi$num.rect.tx[at] <- length(unique(rect.tx))
+  rect.cases <- c(idsRGC_tx, idsRGC_prep_tx, idsRCT_tx, idsRCT_prep_tx)
+  dat$epi$num.rect.cases[at] <- length(unique(rect.cases))
+
+  rect.tx.prep <- c(txRGC_prep, txRCT_prep)
+  dat$epi$num.rect.tx.prep[at] <- length(unique(rect.tx.prep))
+  rect.cases.prep <- c(idsRGC_prep_tx, idsRCT_prep_tx)
+  dat$epi$num.rect.cases.prep[at] <- length(unique(rect.cases.prep))
 
   return(dat)
 }

From eb198fdcb4ddc372a746b8dde069bd4f2ac1860a Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Tue, 27 Sep 2016 20:51:18 -0400
Subject: [PATCH 29/40] Move aging modules back to R from C++

---
 DESCRIPTION                      |  4 +-
 NAMESPACE                        |  1 -
 R/EpiModelHIV-package.R          |  1 -
 R/{RcppExports.R => mod.aging.R} | 39 +++++++++++++++----
 man/aging_het.Rd                 |  4 +-
 man/aging_msm.Rd                 |  4 +-
 src/RcppExports.cpp              | 31 ---------------
 src/mod.aging.cpp                | 66 --------------------------------
 8 files changed, 37 insertions(+), 113 deletions(-)
 rename R/{RcppExports.R => mod.aging.R} (56%)
 delete mode 100644 src/RcppExports.cpp
 delete mode 100644 src/mod.aging.cpp

diff --git a/DESCRIPTION b/DESCRIPTION
index 3aab0fe7..480439f1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -22,13 +22,11 @@ Imports:
     bindata,
     network,
     networkDynamic,
-    Rcpp,
     dplyr
 Suggests:
     knitr,
     testthat
 VignetteBuilder: knitr
-LinkingTo: ergm,
-    Rcpp
+LinkingTo: ergm
 RoxygenNote: 5.0.1
 LazyData: true
diff --git a/NAMESPACE b/NAMESPACE
index 917b82da..ab37e4e7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -70,7 +70,6 @@ import(network)
 import(networkDynamic)
 import(tergm)
 import(tergmLite)
-importFrom(Rcpp,sourceCpp)
 importFrom(dplyr,group_by)
 importFrom(dplyr,summarise)
 importFrom(stats,plogis)
diff --git a/R/EpiModelHIV-package.R b/R/EpiModelHIV-package.R
index 17373f76..dcb8aaee 100644
--- a/R/EpiModelHIV-package.R
+++ b/R/EpiModelHIV-package.R
@@ -20,7 +20,6 @@
 #'
 #' @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
diff --git a/R/RcppExports.R b/R/mod.aging.R
similarity index 56%
rename from R/RcppExports.R
rename to R/mod.aging.R
index 29b12d23..58d949fe 100644
--- a/R/RcppExports.R
+++ b/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)
+}
diff --git a/man/aging_het.Rd b/man/aging_het.Rd
index efbee556..30b7bab7 100644
--- a/man/aging_het.Rd
+++ b/man/aging_het.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/RcppExports.R
+% Please edit documentation in R/mod.aging.R
 \name{aging_het}
 \alias{aging_het}
 \title{Aging Module}
@@ -13,7 +13,7 @@ individual-level attributes, and summary statistics.}
 \item{at}{Current time step.}
 }
 \description{
-This module ages all nodes in the population by one time
+This module ages all active nodes in the population by one time
              unit at each time step.
 }
 \keyword{het}
diff --git a/man/aging_msm.Rd b/man/aging_msm.Rd
index f8b7dbd0..66acb974 100644
--- a/man/aging_msm.Rd
+++ b/man/aging_msm.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/RcppExports.R
+% Please edit documentation in R/mod.aging.R
 \name{aging_msm}
 \alias{aging_msm}
 \title{Aging Module}
@@ -18,7 +18,7 @@ This function returns \code{dat} after updating the nodal attribute
 updated on the three networks.
 }
 \description{
-Module for aging over time for nodes in the population.
+Module for aging over time for active nodes in the population.
 }
 \keyword{module}
 \keyword{msm}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
deleted file mode 100644
index b5ca1d1a..00000000
--- a/src/RcppExports.cpp
+++ /dev/null
@@ -1,31 +0,0 @@
-// Generated by using Rcpp::compileAttributes() -> do not edit by hand
-// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
-
-#include <Rcpp.h>
-
-using namespace Rcpp;
-
-// aging_msm
-List aging_msm(List dat, int at);
-RcppExport SEXP EpiModelHIV_aging_msm(SEXP datSEXP, SEXP atSEXP) {
-BEGIN_RCPP
-    Rcpp::RObject rcpp_result_gen;
-    Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< List >::type dat(datSEXP);
-    Rcpp::traits::input_parameter< int >::type at(atSEXP);
-    rcpp_result_gen = Rcpp::wrap(aging_msm(dat, at));
-    return rcpp_result_gen;
-END_RCPP
-}
-// aging_het
-List aging_het(List dat, int at);
-RcppExport SEXP EpiModelHIV_aging_het(SEXP datSEXP, SEXP atSEXP) {
-BEGIN_RCPP
-    Rcpp::RObject rcpp_result_gen;
-    Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< List >::type dat(datSEXP);
-    Rcpp::traits::input_parameter< int >::type at(atSEXP);
-    rcpp_result_gen = Rcpp::wrap(aging_het(dat, at));
-    return rcpp_result_gen;
-END_RCPP
-}
diff --git a/src/mod.aging.cpp b/src/mod.aging.cpp
deleted file mode 100644
index 013ecc58..00000000
--- a/src/mod.aging.cpp
+++ /dev/null
@@ -1,66 +0,0 @@
-#include <Rcpp.h>
-using namespace Rcpp;
-
-//' @title Aging Module
-//'
-//' @description Module for aging over time for nodes in the population.
-//'
-//' @param dat Master data list object of class \code{dat} containing networks,
-//'        individual-level attributes, and summary statistics.
-//' @param at Current time step.
-//'
-//' @return
-//' This function returns \code{dat} after updating the nodal attribute
-//' \code{age} and \code{sqrt.age}. The \code{sqrt.age} vertex attribute is also
-//' updated on the three networks.
-//'
-//' @keywords module msm
-//' @export
-// [[Rcpp::export]]
-List aging_msm(List dat, int at) {
-
-  List attr = dat["attr"];
-  List param = dat["param"];
-
-  DoubleVector age = as<DoubleVector>(attr["age"]);
-  double tUnit = as<int>(param["time.unit"]);
-
-  DoubleVector nage = age + (tUnit/365);
-
-  attr["age"] = nage;
-  attr["sqrt.age"] = sqrt(nage);
-
-  dat["attr"] = attr;
-
-  return dat;
-}
-
-
-//' @title Aging Module
-//'
-//' @description This module ages all nodes in the population by one time
-//'              unit at each time step.
-//'
-//' @param dat Master data list object of class \code{dat} containing networks,
-//'        individual-level attributes, and summary statistics.
-//' @param at Current time step.
-//'
-//' @keywords module het
-//' @export
-// [[Rcpp::export]]
-List aging_het(List dat, int at) {
-
-  List attr = dat["attr"];
-  List param = dat["param"];
-
-  DoubleVector age = as<DoubleVector>(attr["age"]);
-  double tUnit = as<int>(param["time.unit"]);
-
-  DoubleVector nage = age + (tUnit/365);
-
-  attr["age"] = nage;
-  dat["attr"] = attr;
-
-  return dat;
-}
-

From 729b280254fc80c4bc2a85ec7ea09fa93fc330f2 Mon Sep 17 00:00:00 2001
From: kweiss2 <kevin.weiss@emory.edu>
Date: Fri, 21 Oct 2016 14:42:12 -0400
Subject: [PATCH 30/40] Change parameters and names of parameters to match
 stiPrEP analysis

---
 R/mod.sti.R      | 48 ++++++++++++++++++++++++------------------------
 R/params.R       | 40 ++++++++++++++++++++--------------------
 man/param_msm.Rd | 34 +++++++++++++++++-----------------
 3 files changed, 61 insertions(+), 61 deletions(-)

diff --git a/R/mod.sti.R b/R/mod.sti.R
index f8154b54..19d9b8b9 100644
--- a/R/mod.sti.R
+++ b/R/mod.sti.R
@@ -302,15 +302,15 @@ sti_recov <- function(dat, at) {
 
   # Parameters ----------------------------------------------------------
 
-  rgc.dur.asympt <- dat$param$rgc.dur.asympt
-  ugc.dur.asympt <- dat$param$ugc.dur.asympt
-  gc.dur.tx <- dat$param$gc.dur.tx
-  gc.dur.ntx <- dat$param$gc.dur.ntx
+  rgc.asympt.int <- dat$param$rgc.asympt.int
+  ugc.asympt.int <- dat$param$ugc.asympt.int
+  gc.tx.int <- dat$param$gc.tx.int
+  gc.ntx.int <- dat$param$gc.ntx.int
 
-  rct.dur.asympt <- dat$param$rct.dur.asympt
-  uct.dur.asympt <- dat$param$uct.dur.asympt
-  ct.dur.tx <- dat$param$ct.dur.tx
-  ct.dur.ntx <- dat$param$ct.dur.ntx
+  rct.asympt.int <- dat$param$rct.asympt.int
+  uct.asympt.int <- dat$param$uct.asympt.int
+  ct.tx.int <- dat$param$ct.tx.int
+  ct.ntx.int <- dat$param$ct.ntx.int
 
 
   # GC Recovery ---------------------------------------------------------
@@ -328,9 +328,9 @@ sti_recov <- function(dat, at) {
                              (is.na(dat$attr$uGC.tx.prep) | dat$attr$uGC.tx.prep == 0))
 
   recovRGC_asympt_ntx <- idsRGC_asympt_ntx[which(rbinom(length(idsRGC_asympt_ntx), 1,
-                                                        1/rgc.dur.asympt) == 1)]
+                                                        1/rgc.asympt.int) == 1)]
   recovUGC_asympt_ntx <- idsUGC_asympt_ntx[which(rbinom(length(idsUGC_asympt_ntx), 1,
-                                                        1/ugc.dur.asympt) == 1)]
+                                                        1/ugc.asympt.int) == 1)]
 
   # Symptomatic untreated
   idsRGC_sympt_ntx <- which(dat$attr$rGC == 1 &
@@ -347,14 +347,14 @@ sti_recov <- function(dat, at) {
   # If parameter is null, uses recovery rate of asytomatic untreated
   if (!is.null(gc.dur.ntx)) {
     recovRGC_sympt_ntx <- idsRGC_sympt_ntx[which(rbinom(length(idsRGC_sympt_ntx), 1,
-                                                        1/gc.dur.ntx) == 1)]
+                                                        1/gc.ntx.int) == 1)]
     recovUGC_sympt_ntx <- idsUGC_sympt_ntx[which(rbinom(length(idsUGC_sympt_ntx), 1,
-                                                        1/gc.dur.ntx) == 1)]
+                                                        1/gc.ntx.int) == 1)]
   } else {
     recovRGC_sympt_ntx <- idsRGC_sympt_ntx[which(rbinom(length(idsRGC_sympt_ntx), 1,
-                                                        1/rgc.dur.asympt) == 1)]
+                                                        1/rgc.asympt.int) == 1)]
     recovUGC_sympt_ntx <- idsUGC_sympt_ntx[which(rbinom(length(idsUGC_sympt_ntx), 1,
-                                                        1/ugc.dur.asympt) == 1)]
+                                                        1/ugc.asympt.int) == 1)]
   }
 
   # Treated (asymptomatic and symptomatic)
@@ -366,9 +366,9 @@ sti_recov <- function(dat, at) {
                      (dat$attr$uGC.tx == 1 | dat$attr$uGC.tx.prep == 1))
 
   recovRGC_tx <- idsRGC_tx[which(rbinom(length(idsRGC_tx), 1,
-                                        1/gc.dur.tx) == 1)]
+                                        1/gc.tx.int) == 1)]
   recovUGC_tx <- idsUGC_tx[which(rbinom(length(idsUGC_tx), 1,
-                                        1/gc.dur.tx) == 1)]
+                                        1/gc.tx.int) == 1)]
 
   recovRGC <- c(recovRGC_asympt_ntx, recovRGC_sympt_ntx, recovRGC_tx)
   recovUGC <- c(recovUGC_asympt_ntx, recovUGC_sympt_ntx, recovUGC_tx)
@@ -404,9 +404,9 @@ sti_recov <- function(dat, at) {
                              (is.na(dat$attr$uCT.tx.prep) | dat$attr$uCT.tx.prep == 0))
 
   recovRCT_asympt_ntx <- idsRCT_asympt_ntx[which(rbinom(length(idsRCT_asympt_ntx),
-                                                        1, 1/rct.dur.asympt) == 1)]
+                                                        1, 1/rct.asympt.int) == 1)]
   recovUCT_asympt_ntx <- idsUCT_asympt_ntx[which(rbinom(length(idsUCT_asympt_ntx),
-                                                        1, 1/uct.dur.asympt) == 1)]
+                                                        1, 1/uct.asympt.int) == 1)]
 
   # Symptomatic untreated
   idsRCT_sympt_ntx <- which(dat$attr$rCT == 1 &
@@ -422,14 +422,14 @@ sti_recov <- function(dat, at) {
 
   if (!is.null(ct.dur.ntx)) {
     recovRCT_sympt_ntx <- idsRCT_sympt_ntx[which(rbinom(length(idsRCT_sympt_ntx),
-                                                        1, 1/ct.dur.ntx) == 1)]
+                                                        1, 1/ct.ntx.int) == 1)]
     recovUCT_sympt_ntx <- idsUCT_sympt_ntx[which(rbinom(length(idsUCT_sympt_ntx),
-                                                        1, 1/ct.dur.ntx) == 1)]
+                                                        1, 1/ct.ntx.int) == 1)]
   } else {
     recovRCT_sympt_ntx <- idsRCT_sympt_ntx[which(rbinom(length(idsRCT_sympt_ntx),
-                                                        1, 1/rct.dur.asympt) == 1)]
+                                                        1, 1/rct.asympt.int) == 1)]
     recovUCT_sympt_ntx <- idsUCT_sympt_ntx[which(rbinom(length(idsUCT_sympt_ntx),
-                                                        1, 1/uct.dur.asympt) == 1)]
+                                                        1, 1/uct.asympt.int) == 1)]
   }
 
   # Treated (asymptomatic and symptomatic)
@@ -441,9 +441,9 @@ sti_recov <- function(dat, at) {
                      (dat$attr$uCT.tx == 1 | dat$attr$uCT.tx.prep == 1))
 
   recovRCT_tx <- idsRCT_tx[which(rbinom(length(idsRCT_tx),
-                                        1, 1/ct.dur.tx) == 1)]
+                                        1, 1/ct.tx.int) == 1)]
   recovUCT_tx <- idsUCT_tx[which(rbinom(length(idsUCT_tx),
-                                        1, 1/ct.dur.tx) == 1)]
+                                        1, 1/ct.tx.int) == 1)]
 
 
   recovRCT <- c(recovRCT_asympt_ntx, recovRCT_sympt_ntx, recovRCT_tx)
diff --git a/R/params.R b/R/params.R
index 38bcbe04..e23db1e8 100644
--- a/R/params.R
+++ b/R/params.R
@@ -253,15 +253,15 @@
 #'        chlamydia.
 #' @param uct.sympt.prob Probability of symptoms given infection with urethral
 #'        chlamydia.
-#' @param rgc.dur.asympt Average duration in weeks of asymptomatic rectal gonorrhea.
-#' @param ugc.dur.asympt Average duration in weeks of asymptomatic urethral gonorrhea.
-#' @param gc.dur.tx Average duration in weeks of treated gonorrhea (both sites).
-#' @param gc.dur.ntx Average duration in weeks of untreated, symptomatic gonorrhea (both sites).
+#' @param rgc.asympt.int Average duration in weeks of asymptomatic rectal gonorrhea.
+#' @param ugc.asympt.int Average duration in weeks of asymptomatic urethral gonorrhea.
+#' @param gc.tx.int Average duration in weeks of treated gonorrhea (both sites).
+#' @param gc.ntx.int Average duration in weeks of untreated, symptomatic gonorrhea (both sites).
 #'        If \code{NULL}, uses site-specific durations for asymptomatic infections.
-#' @param rct.dur.asympt Average in weeks duration of asymptomatic rectal chlamydia.
-#' @param uct.dur.asympt Average in weeks duration of asymptomatic urethral chlamydia.
-#' @param ct.dur.tx Average in weeks duration of treated chlamydia (both sites).
-#' @param ct.dur.ntx Average in weeks duration of untreated, symptomatic chlamydia (both sites).
+#' @param rct.asympt.int Average in weeks duration of asymptomatic rectal chlamydia.
+#' @param uct.asympt.int Average in weeks duration of asymptomatic urethral chlamydia.
+#' @param ct.tx.int Average in weeks duration of treated chlamydia (both sites).
+#' @param ct.ntx.int Average in weeks duration of untreated, symptomatic chlamydia (both sites).
 #'        If \code{NULL}, uses site-specific durations for asymptomatic infections.
 #' @param gc.prob.cease Probability of ceasing sexual activity during symptomatic
 #'        infection with gonorrhea.
@@ -425,15 +425,15 @@ param_msm <- function(nwstats,
                       rct.sympt.prob = 0.103517,
                       uct.sympt.prob = 0.885045,
 
-                      rgc.dur.asympt = 34.93723,
-                      ugc.dur.asympt = 36.48328,
-                      gc.dur.tx = 2,
-                      gc.dur.ntx = NULL,
+                      rgc.asympt.int = 35.11851 * 7,
+                      ugc.asympt.int = 35.11851 * 7,
+                      gc.tx.int = 2 * 7,
+                      gc.ntx.int = NULL,
 
-                      rct.dur.asympt = 45.02343,
-                      uct.dur.asympt = 44.88562,
-                      ct.dur.tx = 2,
-                      ct.dur.ntx = NULL,
+                      rct.asympt.int = 44.24538 * 7,
+                      uct.asympt.int = 44.24538 * 7,
+                      ct.tx.int = 2 * 7,
+                      ct.ntx.int = NULL,
 
                       gc.prob.cease = 0,
                       ct.prob.cease = 0,
@@ -449,10 +449,10 @@ param_msm <- function(nwstats,
 
                       sti.cond.rr = 0.3,
 
-                      hiv.rgc.rr = 2.644584,
-                      hiv.ugc.rr = 1.69434,
-                      hiv.rct.rr = 2.644584,
-                      hiv.uct.rr = 1.69434,
+                      hiv.rgc.rr = 2.780673,
+                      hiv.ugc.rr = 1.732363,
+                      hiv.rct.rr = 2.780673,
+                      hiv.uct.rr = 1.732363,
                       hiv.dual.rr = 0.2,
                       ...) {
 
diff --git a/man/param_msm.Rd b/man/param_msm.Rd
index 2c8d8050..84a91dd7 100644
--- a/man/param_msm.Rd
+++ b/man/param_msm.Rd
@@ -52,15 +52,15 @@ param_msm(nwstats, race.method = 1, last.neg.test.B.int = 301,
   rcomp.discl.only = FALSE, rgc.tprob = 0.357698, ugc.tprob = 0.248095,
   rct.tprob = 0.321597, uct.tprob = 0.212965, rgc.sympt.prob = 0.076975,
   ugc.sympt.prob = 0.824368, rct.sympt.prob = 0.103517,
-  uct.sympt.prob = 0.885045, rgc.dur.asympt = 34.93723,
-  ugc.dur.asympt = 36.48328, gc.dur.tx = 2, gc.dur.ntx = NULL,
-  rct.dur.asympt = 45.02343, uct.dur.asympt = 44.88562, ct.dur.tx = 2,
-  ct.dur.ntx = NULL, gc.prob.cease = 0, ct.prob.cease = 0,
-  gc.sympt.prob.tx = 0.9, ct.sympt.prob.tx = 0.85, gc.asympt.prob.tx = 0,
-  ct.asympt.prob.tx = 0, prep.sti.screen.int = 182, prep.sti.prob.tx = 1,
-  prep.continue.stand.tx = TRUE, sti.cond.rr = 0.3, hiv.rgc.rr = 2.644584,
-  hiv.ugc.rr = 1.69434, hiv.rct.rr = 2.644584, hiv.uct.rr = 1.69434,
-  hiv.dual.rr = 0.2, ...)
+  uct.sympt.prob = 0.885045, rgc.asympt.int = 35.11851 * 7,
+  ugc.asympt.int = 35.11851 * 7, gc.tx.int = 2 * 7, gc.ntx.int = NULL,
+  rct.asympt.int = 44.24538 * 7, uct.asympt.int = 44.24538 * 7,
+  ct.tx.int = 2 * 7, ct.ntx.int = NULL, gc.prob.cease = 0,
+  ct.prob.cease = 0, gc.sympt.prob.tx = 0.9, ct.sympt.prob.tx = 0.85,
+  gc.asympt.prob.tx = 0, ct.asympt.prob.tx = 0, prep.sti.screen.int = 182,
+  prep.sti.prob.tx = 1, prep.continue.stand.tx = TRUE, sti.cond.rr = 0.3,
+  hiv.rgc.rr = 2.780673, hiv.ugc.rr = 1.732363, hiv.rct.rr = 2.780673,
+  hiv.uct.rr = 1.732363, hiv.dual.rr = 0.2, ...)
 }
 \arguments{
 \item{nwstats}{Target statistics for the network model. An object of class
@@ -421,22 +421,22 @@ chlamydia.}
 \item{uct.sympt.prob}{Probability of symptoms given infection with urethral
 chlamydia.}
 
-\item{rgc.dur.asympt}{Average duration in weeks of asymptomatic rectal gonorrhea.}
+\item{rgc.asympt.int}{Average duration in weeks of asymptomatic rectal gonorrhea.}
 
-\item{ugc.dur.asympt}{Average duration in weeks of asymptomatic urethral gonorrhea.}
+\item{ugc.asympt.int}{Average duration in weeks of asymptomatic urethral gonorrhea.}
 
-\item{gc.dur.tx}{Average duration in weeks of treated gonorrhea (both sites).}
+\item{gc.tx.int}{Average duration in weeks of treated gonorrhea (both sites).}
 
-\item{gc.dur.ntx}{Average duration in weeks of untreated, symptomatic gonorrhea (both sites).
+\item{gc.ntx.int}{Average duration in weeks of untreated, symptomatic gonorrhea (both sites).
 If \code{NULL}, uses site-specific durations for asymptomatic infections.}
 
-\item{rct.dur.asympt}{Average in weeks duration of asymptomatic rectal chlamydia.}
+\item{rct.asympt.int}{Average in weeks duration of asymptomatic rectal chlamydia.}
 
-\item{uct.dur.asympt}{Average in weeks duration of asymptomatic urethral chlamydia.}
+\item{uct.asympt.int}{Average in weeks duration of asymptomatic urethral chlamydia.}
 
-\item{ct.dur.tx}{Average in weeks duration of treated chlamydia (both sites).}
+\item{ct.tx.int}{Average in weeks duration of treated chlamydia (both sites).}
 
-\item{ct.dur.ntx}{Average in weeks duration of untreated, symptomatic chlamydia (both sites).
+\item{ct.ntx.int}{Average in weeks duration of untreated, symptomatic chlamydia (both sites).
 If \code{NULL}, uses site-specific durations for asymptomatic infections.}
 
 \item{gc.prob.cease}{Probability of ceasing sexual activity during symptomatic

From 3ef3cdb645e955c49c4a3b2acbc253ee45427ab7 Mon Sep 17 00:00:00 2001
From: kweiss2 <kevin.weiss@emory.edu>
Date: Thu, 27 Oct 2016 14:29:48 -0400
Subject: [PATCH 31/40] Fixes to package test (NA for gc.ntx.int, not null)

---
 R/mod.sti.R | 4 ++--
 R/params.R  | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/R/mod.sti.R b/R/mod.sti.R
index 19d9b8b9..21e4cd61 100644
--- a/R/mod.sti.R
+++ b/R/mod.sti.R
@@ -345,7 +345,7 @@ sti_recov <- function(dat, at) {
                             (is.na(dat$attr$uGC.tx.prep) | dat$attr$uGC.tx.prep == 0))
 
   # If parameter is null, uses recovery rate of asytomatic untreated
-  if (!is.null(gc.dur.ntx)) {
+  if (!is.na(gc.ntx.int)) {
     recovRGC_sympt_ntx <- idsRGC_sympt_ntx[which(rbinom(length(idsRGC_sympt_ntx), 1,
                                                         1/gc.ntx.int) == 1)]
     recovUGC_sympt_ntx <- idsUGC_sympt_ntx[which(rbinom(length(idsUGC_sympt_ntx), 1,
@@ -420,7 +420,7 @@ sti_recov <- function(dat, at) {
                             (is.na(dat$attr$uCT.tx) | dat$attr$uCT.tx == 0) &
                             (is.na(dat$attr$uCT.tx.prep) | dat$attr$uCT.tx.prep == 0))
 
-  if (!is.null(ct.dur.ntx)) {
+  if (!is.na(ct.ntx.int)) {
     recovRCT_sympt_ntx <- idsRCT_sympt_ntx[which(rbinom(length(idsRCT_sympt_ntx),
                                                         1, 1/ct.ntx.int) == 1)]
     recovUCT_sympt_ntx <- idsUCT_sympt_ntx[which(rbinom(length(idsUCT_sympt_ntx),
diff --git a/R/params.R b/R/params.R
index e23db1e8..f1468cc9 100644
--- a/R/params.R
+++ b/R/params.R
@@ -428,12 +428,12 @@ param_msm <- function(nwstats,
                       rgc.asympt.int = 35.11851 * 7,
                       ugc.asympt.int = 35.11851 * 7,
                       gc.tx.int = 2 * 7,
-                      gc.ntx.int = NULL,
+                      gc.ntx.int = NA,
 
                       rct.asympt.int = 44.24538 * 7,
                       uct.asympt.int = 44.24538 * 7,
                       ct.tx.int = 2 * 7,
-                      ct.ntx.int = NULL,
+                      ct.ntx.int = NA,
 
                       gc.prob.cease = 0,
                       ct.prob.cease = 0,

From 8da578022bfe2ebb729d3342f74afb40383d0cfa Mon Sep 17 00:00:00 2001
From: kweiss2 <kevin.weiss@emory.edu>
Date: Wed, 2 Nov 2016 10:44:45 -0400
Subject: [PATCH 32/40] Change weeks to days in params.Rd, fix Travis CI issues

---
 R/params.R       | 16 ++++++++--------
 man/param_msm.Rd | 20 ++++++++++----------
 2 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/R/params.R b/R/params.R
index f1468cc9..b7e35a1d 100644
--- a/R/params.R
+++ b/R/params.R
@@ -253,15 +253,15 @@
 #'        chlamydia.
 #' @param uct.sympt.prob Probability of symptoms given infection with urethral
 #'        chlamydia.
-#' @param rgc.asympt.int Average duration in weeks of asymptomatic rectal gonorrhea.
-#' @param ugc.asympt.int Average duration in weeks of asymptomatic urethral gonorrhea.
-#' @param gc.tx.int Average duration in weeks of treated gonorrhea (both sites).
-#' @param gc.ntx.int Average duration in weeks of untreated, symptomatic gonorrhea (both sites).
+#' @param rgc.asympt.int Average duration in days of asymptomatic rectal gonorrhea.
+#' @param ugc.asympt.int Average duration in days of asymptomatic urethral gonorrhea.
+#' @param gc.tx.int Average duration in days of treated gonorrhea (both sites).
+#' @param gc.ntx.int Average duration in days of untreated, symptomatic gonorrhea (both sites).
 #'        If \code{NULL}, uses site-specific durations for asymptomatic infections.
-#' @param rct.asympt.int Average in weeks duration of asymptomatic rectal chlamydia.
-#' @param uct.asympt.int Average in weeks duration of asymptomatic urethral chlamydia.
-#' @param ct.tx.int Average in weeks duration of treated chlamydia (both sites).
-#' @param ct.ntx.int Average in weeks duration of untreated, symptomatic chlamydia (both sites).
+#' @param rct.asympt.int Average in days duration of asymptomatic rectal chlamydia.
+#' @param uct.asympt.int Average in days duration of asymptomatic urethral chlamydia.
+#' @param ct.tx.int Average in days duration of treated chlamydia (both sites).
+#' @param ct.ntx.int Average in days duration of untreated, symptomatic chlamydia (both sites).
 #'        If \code{NULL}, uses site-specific durations for asymptomatic infections.
 #' @param gc.prob.cease Probability of ceasing sexual activity during symptomatic
 #'        infection with gonorrhea.
diff --git a/man/param_msm.Rd b/man/param_msm.Rd
index 84a91dd7..fb71c1fd 100644
--- a/man/param_msm.Rd
+++ b/man/param_msm.Rd
@@ -53,9 +53,9 @@ param_msm(nwstats, race.method = 1, last.neg.test.B.int = 301,
   rct.tprob = 0.321597, uct.tprob = 0.212965, rgc.sympt.prob = 0.076975,
   ugc.sympt.prob = 0.824368, rct.sympt.prob = 0.103517,
   uct.sympt.prob = 0.885045, rgc.asympt.int = 35.11851 * 7,
-  ugc.asympt.int = 35.11851 * 7, gc.tx.int = 2 * 7, gc.ntx.int = NULL,
+  ugc.asympt.int = 35.11851 * 7, gc.tx.int = 2 * 7, gc.ntx.int = NA,
   rct.asympt.int = 44.24538 * 7, uct.asympt.int = 44.24538 * 7,
-  ct.tx.int = 2 * 7, ct.ntx.int = NULL, gc.prob.cease = 0,
+  ct.tx.int = 2 * 7, ct.ntx.int = NA, gc.prob.cease = 0,
   ct.prob.cease = 0, gc.sympt.prob.tx = 0.9, ct.sympt.prob.tx = 0.85,
   gc.asympt.prob.tx = 0, ct.asympt.prob.tx = 0, prep.sti.screen.int = 182,
   prep.sti.prob.tx = 1, prep.continue.stand.tx = TRUE, sti.cond.rr = 0.3,
@@ -421,22 +421,22 @@ chlamydia.}
 \item{uct.sympt.prob}{Probability of symptoms given infection with urethral
 chlamydia.}
 
-\item{rgc.asympt.int}{Average duration in weeks of asymptomatic rectal gonorrhea.}
+\item{rgc.asympt.int}{Average duration in days of asymptomatic rectal gonorrhea.}
 
-\item{ugc.asympt.int}{Average duration in weeks of asymptomatic urethral gonorrhea.}
+\item{ugc.asympt.int}{Average duration in days of asymptomatic urethral gonorrhea.}
 
-\item{gc.tx.int}{Average duration in weeks of treated gonorrhea (both sites).}
+\item{gc.tx.int}{Average duration in days of treated gonorrhea (both sites).}
 
-\item{gc.ntx.int}{Average duration in weeks of untreated, symptomatic gonorrhea (both sites).
+\item{gc.ntx.int}{Average duration in days of untreated, symptomatic gonorrhea (both sites).
 If \code{NULL}, uses site-specific durations for asymptomatic infections.}
 
-\item{rct.asympt.int}{Average in weeks duration of asymptomatic rectal chlamydia.}
+\item{rct.asympt.int}{Average in days duration of asymptomatic rectal chlamydia.}
 
-\item{uct.asympt.int}{Average in weeks duration of asymptomatic urethral chlamydia.}
+\item{uct.asympt.int}{Average in days duration of asymptomatic urethral chlamydia.}
 
-\item{ct.tx.int}{Average in weeks duration of treated chlamydia (both sites).}
+\item{ct.tx.int}{Average in days duration of treated chlamydia (both sites).}
 
-\item{ct.ntx.int}{Average in weeks duration of untreated, symptomatic chlamydia (both sites).
+\item{ct.ntx.int}{Average in days duration of untreated, symptomatic chlamydia (both sites).
 If \code{NULL}, uses site-specific durations for asymptomatic infections.}
 
 \item{gc.prob.cease}{Probability of ceasing sexual activity during symptomatic

From 7ae444409bdcb7c05e071177fc90c9c50b3da106 Mon Sep 17 00:00:00 2001
From: kweiss2 <kevin.weiss@emory.edu>
Date: Wed, 2 Nov 2016 10:47:03 -0400
Subject: [PATCH 33/40] Fix mismatch causing Travis CI to fail

---
 R/params.R       | 4 ++--
 man/param_msm.Rd | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/R/params.R b/R/params.R
index b7e35a1d..212f838b 100644
--- a/R/params.R
+++ b/R/params.R
@@ -257,12 +257,12 @@
 #' @param ugc.asympt.int Average duration in days of asymptomatic urethral gonorrhea.
 #' @param gc.tx.int Average duration in days of treated gonorrhea (both sites).
 #' @param gc.ntx.int Average duration in days of untreated, symptomatic gonorrhea (both sites).
-#'        If \code{NULL}, uses site-specific durations for asymptomatic infections.
+#'        If \code{NA}, uses site-specific durations for asymptomatic infections.
 #' @param rct.asympt.int Average in days duration of asymptomatic rectal chlamydia.
 #' @param uct.asympt.int Average in days duration of asymptomatic urethral chlamydia.
 #' @param ct.tx.int Average in days duration of treated chlamydia (both sites).
 #' @param ct.ntx.int Average in days duration of untreated, symptomatic chlamydia (both sites).
-#'        If \code{NULL}, uses site-specific durations for asymptomatic infections.
+#'        If \code{NA}, uses site-specific durations for asymptomatic infections.
 #' @param gc.prob.cease Probability of ceasing sexual activity during symptomatic
 #'        infection with gonorrhea.
 #' @param ct.prob.cease Probability of ceasing sexual activity during symptomatic
diff --git a/man/param_msm.Rd b/man/param_msm.Rd
index fb71c1fd..73930f21 100644
--- a/man/param_msm.Rd
+++ b/man/param_msm.Rd
@@ -428,7 +428,7 @@ chlamydia.}
 \item{gc.tx.int}{Average duration in days of treated gonorrhea (both sites).}
 
 \item{gc.ntx.int}{Average duration in days of untreated, symptomatic gonorrhea (both sites).
-If \code{NULL}, uses site-specific durations for asymptomatic infections.}
+If \code{NA}, uses site-specific durations for asymptomatic infections.}
 
 \item{rct.asympt.int}{Average in days duration of asymptomatic rectal chlamydia.}
 
@@ -437,7 +437,7 @@ If \code{NULL}, uses site-specific durations for asymptomatic infections.}
 \item{ct.tx.int}{Average in days duration of treated chlamydia (both sites).}
 
 \item{ct.ntx.int}{Average in days duration of untreated, symptomatic chlamydia (both sites).
-If \code{NULL}, uses site-specific durations for asymptomatic infections.}
+If \code{NA}, uses site-specific durations for asymptomatic infections.}
 
 \item{gc.prob.cease}{Probability of ceasing sexual activity during symptomatic
 infection with gonorrhea.}

From a85556e263938d16b4f2bdb13a080e9f08ca9483 Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Sun, 15 Jan 2017 14:47:36 -0500
Subject: [PATCH 34/40] Merge in master to sti-prep branch (new tergmLite
 methods)

---
 R/mod.births.R | 4 ++--
 R/mod.death.R  | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/R/mod.births.R b/R/mod.births.R
index deb25ac5..44026576 100644
--- a/R/mod.births.R
+++ b/R/mod.births.R
@@ -55,7 +55,7 @@ births_msm <- function(dat, at){
   # Update Networks
   if (nBirths > 0) {
     for (i in 1:3) {
-      dat$el[[i]] <- add_vertices(dat$el[[i]], nBirths)
+      dat$el[[i]] <- tergmLite::add_vertices(dat$el[[i]], nBirths)
     }
   }
 
@@ -206,7 +206,7 @@ births_het <- function(dat, at) {
   # Update Population Structure
   if (nBirths > 0) {
     dat <- setBirthAttr_het(dat, at, nBirths)
-    dat$el <- add_vertices(dat$el, nBirths)
+    dat$el <- tergmLite::add_vertices(dat$el, nBirths)
   }
 
   if (unique(sapply(dat$attr, length)) != attributes(dat$el)$n) {
diff --git a/R/mod.death.R b/R/mod.death.R
index 5e4b79bd..fe4cd03d 100644
--- a/R/mod.death.R
+++ b/R/mod.death.R
@@ -160,7 +160,7 @@ deaths_het <- function(dat, at) {
 
   ## 4. Update Population Structure ##
   inactive <- which(dat$attr$active == 0)
-  dat$el <- delete_vertices(dat$el, inactive)
+  dat$el <- tergmLite::delete_vertices(dat$el, inactive)
   dat$attr <- deleteAttr(dat$attr, inactive)
 
   if (unique(sapply(dat$attr, length)) != attributes(dat$el)$n) {

From 8050bc2efff9af80a981bf51929be6024d9cd2e8 Mon Sep 17 00:00:00 2001
From: Samuel Jenness <sjenness@uw.edu>
Date: Sun, 15 Jan 2017 14:47:53 -0500
Subject: [PATCH 35/40] Continue merge

---
 NAMESPACE                    |   1 -
 R/mod.simnet.R               | 289 +++--------------------------------
 README.md                    |   8 +-
 inst/het-test-script.R       |  24 +++
 inst/netsim.mods.R           |  35 +++--
 man/updatenwp_msm.Rd         |  21 ---
 tests/testthat/test-aging.R  |   6 -
 tests/testthat/test-netsim.R |  45 ++++++
 8 files changed, 111 insertions(+), 318 deletions(-)
 create mode 100644 inst/het-test-script.R
 delete mode 100644 man/updatenwp_msm.Rd
 delete mode 100644 tests/testthat/test-aging.R
 create mode 100644 tests/testthat/test-netsim.R

diff --git a/NAMESPACE b/NAMESPACE
index ab37e4e7..617b980f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -57,7 +57,6 @@ export(tx_het)
 export(tx_msm)
 export(update_aiclass_msm)
 export(update_roleclass_msm)
-export(updatenwp_msm)
 export(verbose_het)
 export(verbose_msm)
 export(vl_het)
diff --git a/R/mod.simnet.R b/R/mod.simnet.R
index 16119383..3e7db672 100644
--- a/R/mod.simnet.R
+++ b/R/mod.simnet.R
@@ -19,8 +19,13 @@ simnet_msm <- function(dat, at) {
 
   ## Main network
   nwparam.m <- EpiModel::get_nwparam(dat, network = 1)
-  dat <- updatenwp_msm(dat, network = 1)
 
+  if (dat$param$method == 1) {
+    dat$attr$deg.pers <- get_degree(dat$el[[2]])
+  } else {
+    dat$attr$deg.pers <- paste0(dat$attr$race, get_degree(dat$el[[2]]))
+  }
+  dat <- tergmLite::updateModelTermInputs(dat, network = 1)
 
   dat$el[[1]] <- tergmLite::simulate_network(p = dat$p[[1]],
                                              el = dat$el[[1]],
@@ -28,7 +33,6 @@ simnet_msm <- function(dat, at) {
                                              coef.diss = nwparam.m$coef.diss$coef.adj,
                                              save.changes = TRUE)
 
-
   dat$temp$new.edges <- NULL
   if (at == 2) {
     new.edges.m <- matrix(dat$el[[1]], ncol = 2)
@@ -41,7 +45,13 @@ simnet_msm <- function(dat, at) {
 
   ## Casual network
   nwparam.p <- EpiModel::get_nwparam(dat, network = 2)
-  dat <- updatenwp_msm(dat, network = 2)
+
+  if (dat$param$method == 1) {
+    dat$attr$deg.main <- get_degree(dat$el[[1]])
+  } else {
+    dat$attr$deg.main <- paste0(dat$attr$race, get_degree(dat$el[[1]]))
+  }
+  dat <- tergmLite::updateModelTermInputs(dat, network = 2)
 
   dat$el[[2]] <- tergmLite::simulate_network(p = dat$p[[2]],
                                              el = dat$el[[2]],
@@ -61,7 +71,13 @@ simnet_msm <- function(dat, at) {
 
   ## One-off network
   nwparam.i <- EpiModel::get_nwparam(dat, network = 3)
-  dat <- updatenwp_msm(dat, network = 3)
+
+  if (dat$param$method == 1) {
+    dat$attr$deg.pers <- get_degree(dat$el[[2]])
+  } else {
+    dat$attr$deg.pers <- paste0(dat$attr$race, get_degree(dat$el[[2]]))
+  }
+  dat <- tergmLite::updateModelTermInputs(dat, network = 3)
 
   dat$el[[3]] <- tergmLite::simulate_ergm(p = dat$p[[3]],
                                           el = dat$el[[3]],
@@ -75,6 +91,7 @@ simnet_msm <- function(dat, at) {
 }
 
 
+
 calc_resim_nwstats <- function(dat, at) {
 
   for (nw in 1:3) {
@@ -96,6 +113,7 @@ calc_resim_nwstats <- function(dat, at) {
 }
 
 
+
 #' @title Adjustment for the Edges Coefficient with Changing Network Size
 #'
 #' @description Adjusts the edges coefficients in a dynamic network model
@@ -147,269 +165,6 @@ edges_correct_msm <- function(dat, at) {
 }
 
 
-#' @title Update Network Data Structure and Parameters
-#'
-#' @description Updates the internal data structure containing the main data
-#'              passed into the TERGM resimulation algorithm. This step is
-#'              necessary with the new tergmLite approach.
-#'
-#' @param dat Data object created in initialization module.
-#' @param network Integer value for network number
-#'
-#' @keywords module msm
-#'
-#' @export
-#'
-#'
-updatenwp_msm <- function(dat, network) {
-
-  n <- attributes(dat$el[[1]])$n
-  maxdyads <- choose(n, 2)
-
-  p <- dat$p[[network]]
-
-  if (network == 1) {
-
-    mf <- p$model.form
-    md <- p$model.diss
-    mhf <- p$MHproposal.form
-    mhd <- p$MHproposal.diss
-
-    if (!identical(mf$coef.names, c("edges",
-                                    "nodefactor.deg.pers.1",
-                                    "nodefactor.deg.pers.2",
-                                    "absdiff.sqrt.age",
-                                    "nodematch.role.class.I",
-                                    "nodematch.role.class.R"))) {
-      stop("updatenwp_msm will not currently work with this formula, contact SJ")
-    }
-
-    ## Update model.form ##
-
-    # edges
-    mf$terms[[1]]$maxval <- maxdyads
-
-    # nodefactor("deg.pers")
-    dat$attr$deg.pers <- get_degree(dat$el[[2]])
-
-    nodecov <- dat$attr$deg.pers
-    u <- sort(unique(nodecov))
-    u <- u[-1] # remove base values here
-    nodecov <- match(nodecov, u, nomatch = length(u) + 1)
-    ui <- seq(along = u)
-    inputs <- c(ui, nodecov)
-    mf$terms[[2]]$inputs <- c(length(ui), length(mf$terms[[2]]$coef.names),
-                              length(inputs), inputs)
-
-    # absdiff("sqrt.age")
-    nodecov <- dat$attr$sqrt.age
-    power <- 1
-    inputs <- c(power, nodecov)
-    mf$terms[[3]]$inputs <- c(0, length(mf$terms[[3]]$coef.names),
-                              length(inputs), inputs)
-
-    # nodematch("role.class)
-    nodecov <- dat$attr$role.class
-    u <- sort(unique(nodecov))
-    u <- u[1:2] # keep = 1:2
-    nodecov <- match(nodecov, u, nomatch = length(u) + 1)
-    dontmatch <- nodecov == (length(u) + 1)
-    nodecov[dontmatch] <- length(u) + (1:sum(dontmatch))
-    ui <- seq(along = u)
-    inputs <- c(ui, nodecov)
-    mf$terms[[4]]$inputs <- c(0, length(mf$terms[[4]]$coef.names),
-                              length(inputs), inputs)
-
-    ## combined maxval ##
-    mf$maxval[1] <- maxdyads
-
-
-    ## Update model.diss ##
-    md$terms[[1]]$maxval <- maxdyads
-    md$maxval <- maxdyads
-
-    ## Update MHproposal.form ##
-    mhf$arguments$constraints$bd$attribs <-
-           matrix(rep(mhf$arguments$constraints$bd$attribs[1], n), ncol = 1)
-    mhf$arguments$constraints$bd$maxout <-
-            matrix(rep(mhf$arguments$constraints$bd$maxout[1], n), ncol = 1)
-    mhf$arguments$constraints$bd$maxin <- matrix(rep(n - 1, n), ncol = 1)
-    mhf$arguments$constraints$bd$minout <-
-           mhf$arguments$constraints$bd$minin <- matrix(rep(0, n), ncol = 1)
-
-    ## Update MHproposal.diss ##
-    mhd$arguments$constraints$bd <- mhf$arguments$constraints$bd
-
-    dat$p[[network]] <- list(model.form = mf, model.diss = md,
-                             MHproposal.form = mhf, MHproposal.diss = mhd)
-
-  }
-
-  if (network == 2) {
-
-    mf <- p$model.form
-    md <- p$model.diss
-    mhf <- p$MHproposal.form
-    mhd <- p$MHproposal.diss
-
-    if (!identical(mf$coef.names, c("edges",
-                                    "nodefactor.deg.main.1",
-                                    "concurrent",
-                                    "absdiff.sqrt.age",
-                                    "nodematch.role.class.I",
-                                    "nodematch.role.class.R"))) {
-      stop("updatenwp_msm will not currently work with this formula, contact SJ")
-    }
-
-
-    ## Update model.form ##
-
-    # edges
-    mf$terms[[1]]$maxval <- maxdyads
-
-    # nodefactor("deg.main")
-    dat$attr$deg.main <- get_degree(dat$el[[1]])
-
-    nodecov <- dat$attr$deg.main
-    u <- sort(unique(nodecov))
-    u <- u[-1] # remove base values here
-    nodecov <- match(nodecov, u, nomatch = length(u) + 1)
-    ui <- seq(along = u)
-    inputs <- c(ui, nodecov)
-    mf$terms[[2]]$inputs <- c(length(ui), length(mf$terms[[2]]$coef.names),
-                              length(inputs), inputs)
-
-    # concurrent
-    mf$terms[[3]]$maxval <- n
-
-    # absdiff("sqrt.age")
-    nodecov <- dat$attr$sqrt.age
-    power <- 1
-    inputs <- c(power, nodecov)
-    mf$terms[[4]]$inputs <- c(0, length(mf$terms[[4]]$coef.names),
-                              length(inputs), inputs)
-
-    # nodematch("role.class)
-    nodecov <- dat$attr$role.class
-    u <- sort(unique(nodecov))
-    u <- u[1:2] # keep = 1:2
-    nodecov <- match(nodecov, u, nomatch = length(u) + 1)
-    dontmatch <- nodecov == (length(u) + 1)
-    nodecov[dontmatch] <- length(u) + (1:sum(dontmatch))
-    ui <- seq(along = u)
-    inputs <- c(ui, nodecov)
-    mf$terms[[5]]$inputs <- c(0, length(mf$terms[[5]]$coef.names),
-                              length(inputs), inputs)
-
-    ## combined maxval ##
-    mf$maxval[1] <- maxdyads
-    mf$maxval[3] <- n
-
-    ## Update model.diss ##
-    md$terms[[1]]$maxval <- maxdyads
-    md$maxval <- maxdyads
-
-    ## Update MHproposal.form ##
-    mhf$arguments$constraints$bd$attribs <-
-           matrix(rep(mhf$arguments$constraints$bd$attribs[1], n), ncol = 1)
-    mhf$arguments$constraints$bd$maxout <-
-            matrix(rep(mhf$arguments$constraints$bd$maxout[1], n), ncol = 1)
-    mhf$arguments$constraints$bd$maxin <- matrix(rep(n - 1, n), ncol = 1)
-    mhf$arguments$constraints$bd$minout <-
-           mhf$arguments$constraints$bd$minin <- matrix(rep(0, n), ncol = 1)
-
-    ## Update MHproposal.diss ##
-    mhd$arguments$constraints$bd <- mhf$arguments$constraints$bd
-
-    dat$p[[network]] <- list(model.form = mf, model.diss = md,
-                             MHproposal.form = mhf, MHproposal.diss = mhd)
-
-  }
-
-  if (network == 3) {
-
-    mf <- p$model.form
-    mhf <- p$MHproposal
-
-    if (!identical(mf$coef.names, c("edges",
-                                    "nodefactor.deg.main.deg.pers.0.1",
-                                    "nodefactor.deg.main.deg.pers.0.2",
-                                    "nodefactor.deg.main.deg.pers.1.0",
-                                    "nodefactor.deg.main.deg.pers.1.1",
-                                    "nodefactor.deg.main.deg.pers.1.2",
-                                    "nodefactor.riskg.1",
-                                    "nodefactor.riskg.2",
-                                    "nodefactor.riskg.4",
-                                    "nodefactor.riskg.5",
-                                    "absdiff.sqrt.age",
-                                    "nodematch.role.class.I",
-                                    "nodematch.role.class.R"))) {
-      stop("updatenwp_msm will not currently work with this formula, contact SJ")
-    }
-
-
-    ## Update model.form ##
-
-    # edges
-    mf$terms[[1]]$maxval <- maxdyads
-
-    # nodefactor(c("deg.main", "deg.pers"))
-    # current main degree already written in last conditional block
-    dat$attr$deg.pers <- get_degree(dat$el[[2]])
-
-    nodecov <- do.call(paste, c(sapply(c("deg.main", "deg.pers"),
-                                       function(oneattr) dat$attr[[oneattr]],
-                                       simplify = FALSE), sep = "."))
-    u <- sort(unique(nodecov))
-    u <- u[-1] # remove base values here
-    nodecov <- match(nodecov, u, nomatch = length(u) + 1)
-    ui <- seq(along = u)
-    inputs <- c(ui, nodecov)
-    mf$terms[[2]]$inputs <- c(length(ui), length(mf$terms[[2]]$coef.names),
-                              length(inputs), inputs)
-
-
-    # nodefactor("riskg", base = 3)
-    nodecov <- dat$attr$riskg
-    u <- sort(unique(nodecov))
-    u <- u[-3] # remove base values here
-    nodecov <- match(nodecov, u, nomatch = length(u) + 1)
-    ui <- seq(along = u)
-    inputs <- c(ui, nodecov)
-    mf$terms[[3]]$inputs <- c(length(ui), length(mf$terms[[3]]$coef.names),
-                              length(inputs), inputs)
-
-    # absdiff("sqrt.age")
-    nodecov <- dat$attr$sqrt.age
-    power <- 1
-    inputs <- c(power, nodecov)
-    mf$terms[[4]]$inputs <- c(0, length(mf$terms[[4]]$coef.names),
-                              length(inputs), inputs)
-
-    # nodematch("role.class)
-    nodecov <- dat$attr$role.class
-    u <- sort(unique(nodecov))
-    u <- u[1:2] # keep = 1:2
-    nodecov <- match(nodecov, u, nomatch = length(u) + 1)
-    dontmatch <- nodecov == (length(u) + 1)
-    nodecov[dontmatch] <- length(u) + (1:sum(dontmatch))
-    ui <- seq(along = u)
-    inputs <- c(ui, nodecov)
-    mf$terms[[5]]$inputs <- c(0, length(mf$terms[[5]]$coef.names),
-                              length(inputs), inputs)
-
-    ## combined maxval ##
-    mf$maxval[1] <- maxdyads
-
-    ## Update MHproposal ##
-    # no changes
-
-    dat$p[[network]] <- list(model.form = mf, MHproposal = mhf)
-  }
-
-  return(dat)
-}
-
 
 
 # HET -----------------------------------------------------------------
diff --git a/README.md b/README.md
index 8046f6ce..7df42c3e 100644
--- a/README.md
+++ b/README.md
@@ -20,9 +20,6 @@ devtools::install_github("statnet/EpiModelHIV")
 
 Documentation on using this software package is forthcoming, although limited function documentation is provided within the package and available with the `help(package = "EpiModelHIV")` command.
 
-### Note on Repository and Package Name
-This Github repository `EpiModelHIV` and the R package within it were previously named `EpiModelHIVmsm`. On 6/24/2016, we merged that MSM package with our `EpiModelHIVhet` (HIV models for heterosexuals) into this combined repository and package. All scripts from those separate packages should still function with this current version after changing the input to `library`. 
-
 ## Software Development Team
 
 Author                                                        | Department                 | Institution
@@ -35,7 +32,8 @@ Author                                                        | Department
 
 `EpiModelHIV` has been used in the following scientific articles:
 
-1. Jenness SM, Goodreau SM, Morris M, Cassels S. Effectiveness of Combination Packages for HIV-1 Prevention in Sub-Saharan Africa Depends on Partnership Network Structure. _Sexually Transmitted Infections._ [DOI: 10.1136/sextrans-2015-052476.](http://sti.bmj.com/content/early/2016/06/09/sextrans-2015-052476.abstract)
+1. Jenness SM, Goodreau SM, Morris M, Cassels S. Effectiveness of Combination Packages for HIV-1 Prevention in Sub-Saharan Africa Depends on Partnership Network Structure. _Sexually Transmitted Infections._ 2016; 92(8): 619-624. [LINK](http://sti.bmj.com/content/early/2016/06/09/sextrans-2015-052476.abstract)
 
-2. Jenness SM, Goodreau SM, Rosenberg E, Beylerian EN, Hoover KW, Smith DK, Sullivan P. Impact of CDC’s HIV Preexposure Prophylaxis Guidelines among MSM in the United States. _Journal of Infectious Diseases._ [DOI: 10.1093/infdis/jiw223.](http://jid.oxfordjournals.org/content/early/2016/07/12/infdis.jiw223.full)
+2. Jenness SM, Goodreau SM, Rosenberg E, Beylerian EN, Hoover KW, Smith DK, Sullivan P. Impact of CDC’s HIV Preexposure Prophylaxis Guidelines among MSM in the United States. _Journal of Infectious Diseases._ 2016; 214(12): 1800-1807. [LINK](http://jid.oxfordjournals.org/content/early/2016/07/12/infdis.jiw223.full)
 
+3. Jenness SM, Sharma A, Goodreau SM, Rosenberg ES, Weiss KM, Hoover KW, Smith DK, Sullivan P. Individual HIV Risk versus Population Impact of Risk Compensation after HIV Preexposure Prophylaxis Initiation among Men Who Have Sex with Men. _PLoS One._ 2017; 12(1): e0169484. [LINK](http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0169484)
diff --git a/inst/het-test-script.R b/inst/het-test-script.R
new file mode 100644
index 00000000..fa721270
--- /dev/null
+++ b/inst/het-test-script.R
@@ -0,0 +1,24 @@
+
+## Heterosexual model test script
+
+library(EpiModelHIV)
+
+st <- make_nw_het(part.dur = 2013)
+est <- netest(st$nw,
+              formation = st$formation,
+              target.stats = st$stats,
+              coef.form = -Inf,
+              coef.diss = st$coef.diss,
+              constraints = ~bd(maxout = 3),
+              set.control.ergm = control.ergm(MCMLE.maxit = 500, MPLE.type = "penalized"))
+
+dx <- netdx(est, nsims = 5, nsteps = 250,
+            set.control.ergm = control.simulate.ergm(MCMC.burnin = 1e6))
+print(dx)
+plot(dx)
+
+param <- param_het()
+init <- init_het(i.prev.male = 0.25, i.prev.feml = 0.25)
+control <- control_het(nsteps = 2600)
+
+sim <- netsim(est, param, init, control)
diff --git a/inst/netsim.mods.R b/inst/netsim.mods.R
index 75e38a96..9a96a237 100644
--- a/inst/netsim.mods.R
+++ b/inst/netsim.mods.R
@@ -1,31 +1,34 @@
-rm(list=ls())
+rm(list = ls())
 suppressMessages(library("EpiModelHIV"))
-sourceDir("R/")
 
 data(est)
 data(st)
-est
-st
-param <- param_msm(nwstats = st, 
+
+param <- param_msm(nwstats = st,
                    ai.scale = 1.323,
                    prep.coverage = 0)
-init <- init_msm(nwstats = st, 
-                 prev.B = 0.253, 
+init <- init_msm(nwstats = st,
+                 prev.B = 0.253,
                  prev.W = 0.253)
-control <- control_msm(simno = 0.253, 
+control <- control_msm(simno = 0.253,
                        nsteps = 52,
-                       nsims = 5, 
-                       ncores = 1, 
+                       nsims = 5,
+                       ncores = 1,
                        save.nwstats = TRUE,
                        verbose.int = 1)
-sim <- netsim(est, param, init, control)
+# sim <- netsim(est, param, init, control)
 
+# debug(stergm_prep)
 
 
 at <- 1
 dat <- initialize_msm(est, param, init, control, s = 1)
 # dat <- reinit_msm(sim, param, init, control, s = 1)
 
+# mf <- dat$p[[1]]$model.form
+# mf$terms[[4]]
+
+
 at <- at + 1
 dat <- aging_msm(dat, at)
 dat <- deaths_msm(dat, at)
@@ -34,11 +37,9 @@ dat <- test_msm(dat, at)
 dat <- tx_msm(dat, at)
 dat <- prep_msm(dat, at)
 dat <- progress_msm(dat, at)
-dat <- update_vl_msm(dat, at)
-dat <- update_aiclass_msm(dat, at)
-dat <- update_roleclass_msm(dat, at)
-dat <- edges_correct_msm(dat, at)
-dat <- updatenwp_msm(dat, at)
+dat <- vl_msm(dat, at)
+# dat <- update_aiclass_msm(dat, at)
+# dat <- update_roleclass_msm(dat, at)
 dat <- simnet_msm(dat, at)
 dat <- disclose_msm(dat, at)
 dat <- acts_msm(dat, at)
@@ -47,5 +48,3 @@ dat <- riskhist_msm(dat, at)
 dat <- position_msm(dat, at)
 dat <- trans_msm(dat, at)
 dat <- prevalence_msm(dat, at)
-verbose_msm(dat, type = "progress", s = 1, at)
-
diff --git a/man/updatenwp_msm.Rd b/man/updatenwp_msm.Rd
deleted file mode 100644
index 903ff227..00000000
--- a/man/updatenwp_msm.Rd
+++ /dev/null
@@ -1,21 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/mod.simnet.R
-\name{updatenwp_msm}
-\alias{updatenwp_msm}
-\title{Update Network Data Structure and Parameters}
-\usage{
-updatenwp_msm(dat, network)
-}
-\arguments{
-\item{dat}{Data object created in initialization module.}
-
-\item{network}{Integer value for network number}
-}
-\description{
-Updates the internal data structure containing the main data
-             passed into the TERGM resimulation algorithm. This step is
-             necessary with the new tergmLite approach.
-}
-\keyword{module}
-\keyword{msm}
-
diff --git a/tests/testthat/test-aging.R b/tests/testthat/test-aging.R
deleted file mode 100644
index c4155a38..00000000
--- a/tests/testthat/test-aging.R
+++ /dev/null
@@ -1,6 +0,0 @@
-context("Aging Module")
-
-test_that("Aging module functions properly", {
-
-  NULL
-})
\ No newline at end of file
diff --git a/tests/testthat/test-netsim.R b/tests/testthat/test-netsim.R
new file mode 100644
index 00000000..1472c176
--- /dev/null
+++ b/tests/testthat/test-netsim.R
@@ -0,0 +1,45 @@
+context("Model Runs")
+
+test_that("Burnin model", {
+
+  data(st)
+  data(est)
+
+  param <- param_msm(nwstats = st)
+  init <- init_msm(nwstats = st)
+  control <- control_msm(simno = 1, nsteps = 10, verbose = FALSE)
+
+  sim <- netsim(est, param, init, control)
+  expect_is(sim, "netsim")
+
+})
+
+test_that("Follow-up model", {
+
+  data(st)
+  data(est)
+
+  param <- param_msm(nwstats = st,
+                     prep.start = 10)
+  init <- init_msm(nwstats = st)
+  control <- control_msm(simno = 1, nsteps = 10,
+                         save.other = c("attr", "temp", "riskh", "el", "p"),
+                         verbose = FALSE)
+
+  sim <- netsim(est, param, init, control)
+
+  param <- param_msm(nwstats = st,
+                     prep.start = 10,
+                     prep.elig.model = "cdc3",
+                     prep.coverage = 0.5,
+                     prep.risk.int = 182,
+                     prep.class.prob = reallocate_pcp(reall = 0),
+                     prep.class.hr = c(1, 0.69, 0.19, 0.05))
+  init <- init_msm(nwstats = st)
+  control <- control_msm(simno = 1, start = 11, nsteps = 20,
+                         verbose = FALSE, initialize.FUN = reinit_msm)
+
+  sim2 <- netsim(sim, param, init, control)
+  expect_is(sim2, "netsim")
+
+})

From b7ad30420341a6faa70656daadfafd475f8848b0 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Thu, 4 May 2017 12:44:45 -0400
Subject: [PATCH 36/40] Remove RcppExports.cpp

---
 src/RcppExports.cpp | 31 -------------------------------
 1 file changed, 31 deletions(-)
 delete mode 100644 src/RcppExports.cpp

diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
deleted file mode 100644
index 7d2ad9b6..00000000
--- a/src/RcppExports.cpp
+++ /dev/null
@@ -1,31 +0,0 @@
-// This file was generated by Rcpp::compileAttributes
-// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
-
-#include <Rcpp.h>
-
-using namespace Rcpp;
-
-// aging_msm
-List aging_msm(List dat, int at);
-RcppExport SEXP EpiModelHIV_aging_msm(SEXP datSEXP, SEXP atSEXP) {
-BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
-    Rcpp::traits::input_parameter< List >::type dat(datSEXP);
-    Rcpp::traits::input_parameter< int >::type at(atSEXP);
-    __result = Rcpp::wrap(aging_msm(dat, at));
-    return __result;
-END_RCPP
-}
-// aging_het
-List aging_het(List dat, int at);
-RcppExport SEXP EpiModelHIV_aging_het(SEXP datSEXP, SEXP atSEXP) {
-BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
-    Rcpp::traits::input_parameter< List >::type dat(datSEXP);
-    Rcpp::traits::input_parameter< int >::type at(atSEXP);
-    __result = Rcpp::wrap(aging_het(dat, at));
-    return __result;
-END_RCPP
-}

From 6f9f93b95f856ec4164e7994bae31a28ee3953af Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Thu, 4 May 2017 12:47:41 -0400
Subject: [PATCH 37/40] Update documentation

---
 DESCRIPTION             | 9 +++++----
 R/EpiModelHIV-package.R | 6 +++---
 2 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 7f8b067c..e395e1d7 100644
--- a/DESCRIPTION
+++ b/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
diff --git a/R/EpiModelHIV-package.R b/R/EpiModelHIV-package.R
index dcb8aaee..cdbff24d 100644
--- a/R/EpiModelHIV-package.R
+++ b/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
 #' }

From eef93ddc97d5f1b18520c5450e2a59b08517a29e Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Thu, 4 May 2017 12:47:51 -0400
Subject: [PATCH 38/40] Build docs

---
 man/EpiModelHIV-package.Rd | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/man/EpiModelHIV-package.Rd b/man/EpiModelHIV-package.Rd
index 252d145f..ab477c4b 100644
--- a/man/EpiModelHIV-package.Rd
+++ b/man/EpiModelHIV-package.Rd
@@ -4,13 +4,13 @@
 \name{EpiModelHIV-package}
 \alias{EpiModelHIV-package}
 \alias{EpiModelHIV}
-\title{HIV Transmission Dynamics among MSM and Heterosexuals}
+\title{Network-Based Epidemic Modeling of HIV Transmission among MSM and Heterosexual Populations}
 \description{
 \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
 }

From 05c9274abb0b6b1b878a3de1c383e3f4471ed398 Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Thu, 4 May 2017 12:50:41 -0400
Subject: [PATCH 39/40] Update readme

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 0d270312..dfed8dd3 100644
--- a/README.md
+++ b/README.md
@@ -15,6 +15,7 @@ You can install `EpiModelHIV` in R using `devtools`:
 ```
 install.packages("EpiModel", dependencies = TRUE)
 devtools::install_github("statnet/tergmLite")
+devtools::install_github("statnet/EpiModelHPC")
 devtools::install_github("statnet/EpiModelHIV")
 ```
 
@@ -34,7 +35,6 @@ Kevin M. Weiss                                                | Department of Ep
 `EpiModelHIV` has been used in the following scientific articles:
 
 1. Jenness SM, Goodreau SM, Morris M, Cassels S. Effectiveness of Combination Packages for HIV-1 Prevention in Sub-Saharan Africa Depends on Partnership Network Structure. _Sexually Transmitted Infections._ 2016; 92(8): 619-624. [LINK](http://sti.bmj.com/content/early/2016/06/09/sextrans-2015-052476.abstract)
-<<<<<<< HEAD
 
 2. Jenness SM, Goodreau SM, Rosenberg E, Beylerian EN, Hoover KW, Smith DK, Sullivan P. Impact of CDC’s HIV Preexposure Prophylaxis Guidelines among MSM in the United States. _Journal of Infectious Diseases._ 2016; 214(12): 1800-1807. [LINK](http://jid.oxfordjournals.org/content/early/2016/07/12/infdis.jiw223.full)
 

From 7ccebdb55308cf910b761efa192adade38ee136f Mon Sep 17 00:00:00 2001
From: smjenness <sjenness@uw.edu>
Date: Thu, 4 May 2017 12:51:45 -0400
Subject: [PATCH 40/40] Update citations

---
 README.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/README.md b/README.md
index dfed8dd3..a7338810 100644
--- a/README.md
+++ b/README.md
@@ -42,3 +42,4 @@ Kevin M. Weiss                                                | Department of Ep
 
 4. Goodreau SM, Rosenberg ES, Jenness SM, Luisi N, Stansfield SE, Millett G, Sullivan P. Sources of Racial Disparities in HIV Prevalence among Men Who Have Sex with Men in Atlanta, GA: A Modeling Study. _Lancet HIV._ Epub ahead of print. DOI: 10.1016/ S2352-3018(17)30067. [LINK](https://www.ncbi.nlm.nih.gov/pubmed/28431923)
 
+5. Jenness SM, Weiss KM, Goodreau SM, Rosenberg E, Gift T, Chesson H, Hoover KW, Smith DK, Liu AY, Sullivan P. Incidence of Gonorrhea and Chlamydia Following HIV Preexposure Prophylaxis among Men Who Have Sex with Men: A Modeling Study. _Clinical Infectious Diseases._ In Press.