Skip to content

Commit

Permalink
Add branch files
Browse files Browse the repository at this point in the history
  • Loading branch information
@smjenness
smjenness committed Jun 22, 2021
1 parent 26ea27a commit 962e109
Show file tree
Hide file tree
Showing 5 changed files with 168 additions and 83 deletions.
34 changes: 3 additions & 31 deletions R/mod.hivtrans.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,6 @@ hivtrans_msm <- function(dat, at) {
rCT <- dat$attr$rCT
uCT <- dat$attr$uCT
race <- dat$attr$race
tx.status <- dat$attr$tx.status

# Parameters
URAI.prob <- dat$param$URAI.prob
Expand Down Expand Up @@ -89,7 +88,6 @@ hivtrans_msm <- function(dat, at) {
# Attributes of infected
ip.vl <- vl[disc.ip[, 1]]
ip.stage <- stage[disc.ip[, 1]]
ip.txStat <- tx.status[disc.ip[, 1]]

# Attributes of susceptible
ip.prep <- prepStat[disc.ip[, 2]]
Expand All @@ -100,10 +98,6 @@ hivtrans_msm <- function(dat, at) {
# Base TP from VL
ip.tprob <- pmin(0.99, URAI.prob * 2.45^(ip.vl - 4.5))

# Adjustment (based on Supervie JAIDS) for VL Suppressed, on ART
ip.noTrans <- which(ip.vl <= log10(200) & ip.txStat == 1)
ip.tprob[ip.noTrans] <- 2.2/1e5

# Transform to log odds
ip.tlo <- log(ip.tprob/(1 - ip.tprob))

Expand Down Expand Up @@ -151,7 +145,6 @@ hivtrans_msm <- function(dat, at) {
# Attributes of infected
rp.vl <- vl[disc.rp[, 2]]
rp.stage <- stage[disc.rp[, 2]]
rp.txStat <- tx.status[disc.rp[, 2]]

# Attributes of susceptible
rp.circ <- circ[disc.rp[, 1]]
Expand All @@ -163,10 +156,6 @@ hivtrans_msm <- function(dat, at) {
# Base TP from VL
rp.tprob <- pmin(0.99, UIAI.prob * 2.45^(rp.vl - 4.5))

# Adjustment (based on Supervie JAIDS) for VL Suppressed, on ART
rp.noTrans <- which(rp.vl <= log10(200) & rp.txStat == 1)
rp.tprob[rp.noTrans] <- 2.2/1e5

# Transform to log odds
rp.tlo <- log(rp.tprob/(1 - rp.tprob))

Expand Down Expand Up @@ -237,37 +226,20 @@ hivtrans_msm <- function(dat, at) {
dat$attr$cuml.time.off.tx[infected] <- 0

# Attributes of transmitter
transmitter <- as.numeric(c(disc.ip[trans.ip == 1, 1],
disc.rp[trans.rp == 1, 2]))
transmitter <- c(disc.ip[trans.ip == 1, 1],
disc.rp[trans.rp == 1, 2])
tab.trans <- table(transmitter)
uni.trans <- as.numeric(names(tab.trans))
dat$attr$count.trans[uni.trans] <- dat$attr$count.trans[uni.trans] +
as.numeric(tab.trans)
}
}

# Summary Output
dat$epi$incid[at] <- length(infected)
dat$epi$incid.B[at] <- sum(dat$attr$race[infected] == 1)
dat$epi$incid.H[at] <- sum(dat$attr$race[infected] == 2)
dat$epi$incid.W[at] <- sum(dat$attr$race[infected] == 3)

if (length(infected) > 0) {
dat$epi$incid.undx[at] <- sum(dat$attr$diag.status[transmitter] == 0)
dat$epi$incid.dx[at] <- sum(dat$attr$diag.status[transmitter] == 1 &
dat$attr$cuml.time.on.tx[transmitter] == 0)
dat$epi$incid.linked[at] <- sum(dat$attr$diag.status[transmitter] == 1 &
dat$attr$cuml.time.on.tx[transmitter] > 0 &
dat$attr$vl[transmitter] > log10(200))
dat$epi$incid.vsupp[at] <- sum(dat$attr$diag.status[transmitter] == 1 &
dat$attr$cuml.time.on.tx[transmitter] > 0 &
dat$attr$vl[transmitter] <= log10(200))
} else {
dat$epi$incid.undx[at] <- 0
dat$epi$incid.dx[at] <- 0
dat$epi$incid.linked[at] <- 0
dat$epi$incid.vsupp[at] <- 0
}

return(dat)
}

Expand Down
158 changes: 144 additions & 14 deletions R/mod.prep.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,6 @@ prep_msm <- function(dat, at) {

# Core Attributes
active <- dat$attr$active
status <- dat$attr$status
diag.status <- dat$attr$diag.status
lnt <- dat$attr$last.neg.test

Expand All @@ -41,6 +40,16 @@ prep_msm <- function(dat, at) {
prepStartTime <- dat$attr$prepStartTime
prepLastStiScreen <- dat$attr$prepLastStiScreen

if (is.null(dat$attr$pOptimInitStatus)) {
n <- length(active)
dat$attr$pOptimInitStatus <- rep(NA, n)
dat$attr$pOptimAdhrStatus <- rep(NA, n)
dat$attr$pOptimRetnStatus <- rep(NA, n)
}
pOptimInitStatus <- dat$attr$pOptimInitStatus
pOptimAdhrStatus <- dat$attr$pOptimAdhrStatus
pOptimRetnStatus <- dat$attr$pOptimRetnStatus


# Parameters --------------------------------------------------------------

Expand All @@ -50,11 +59,18 @@ prep_msm <- function(dat, at) {
prep.risk.reassess.method <- dat$param$prep.risk.reassess.method
prep.discont.rate <- dat$param$prep.discont.rate

prep.optim.start <- dat$param$prep.optim.start
prep.optim.init.prob <- dat$param$prep.optim.init.prob
prep.start.prob.optim <- dat$param$prep.start.prob.optim
prep.adhr.dist.optim <- dat$param$prep.adhr.dist.optim
prep.optim.adhr.cap <- dat$param$prep.optim.adhr.cap
prep.optim.retn.cap <- dat$param$prep.optim.retn.cap
prep.discont.rate.optim <- dat$param$prep.discont.rate.optim


# Indications -------------------------------------------------------------

ind1 <- dat$attr$prep.ind.uai.mono
# ind2 <- dat$attr$prep.ind.uai.nmain
ind2 <- dat$attr$prep.ind.uai.conc
ind3 <- dat$attr$prep.ind.sti

Expand All @@ -72,6 +88,8 @@ prep_msm <- function(dat, at) {
base.cond.yes <- which(active == 1 & diag.status == 0)
idsIndic <- intersect(idsIndic, base.cond.yes)

idsNewIndic <- intersect(which(prepElig == 0), idsIndic)

# Set eligibility to 1 if indications
prepElig[idsIndic] <- 1

Expand All @@ -98,54 +116,147 @@ prep_msm <- function(dat, at) {
idsStpInd <- intersect(idsNoIndic, idsRiskAssess)
}

# Random discontinuation
idsEligStpRand <- which(active == 1 & prepStat == 1)
# Random Discontinuation (Non-Intervention)
idsEligStpRand <- which(active == 1 & prepStat == 1 & is.na(pOptimRetnStatus))
vecStpRand <- rbinom(length(idsEligStpRand), 1, prep.discont.rate)
idsStpRand <- idsEligStpRand[which(vecStpRand == 1)]

# Random Discontinuation (Intervention)
idsStpRandOptim <- NULL
if (at >= prep.optim.start) {
idsEligStpRandOptim <- which(active == 1 & prepStat == 1 & pOptimRetnStatus == 1)
vecStpRandOptim <- rbinom(length(idsEligStpRandOptim), 1, prep.discont.rate.optim)
idsStpRandOptim <- idsEligStpRandOptim[which(vecStpRandOptim == 1)]
}

# Diagnosis
idsStpDx <- which(active == 1 & prepStat == 1 & diag.status == 1)

# Death
idsStpDth <- which(active == 0 & prepStat == 1)

# Reset PrEP status
idsStp <- c(idsStpInd, idsStpRand, idsStpDx, idsStpDth)
idsStp <- c(idsStpInd, idsStpRand, idsStpRandOptim, idsStpDx, idsStpDth)

# Update attributes for stoppers
prepStat[idsStp] <- 0
prepLastRisk[idsStp] <- NA
prepStartTime[idsStp] <- NA
prepLastStiScreen[idsStp] <- NA
prepClass[idsStp] <- NA

PrEPStopsInd <- length(idsStpInd)
PrEPStopsRand <- length(idsStpRand)
PrEPStopsRandOptim <- length(idsStpRandOptim)

## Initiation ----------------------------------------------------------------

## Eligibility ##
## Initiation ----------------------------------------------------------------

# Indications to start
## Initiation (non-intervention)
if (prep.require.lnt == TRUE) {
idsEligStart <- which(prepStat == 0 & lnt == at)
} else {
idsEligStart <- which(prepStat == 0)
}

idsEligStart <- intersect(idsIndic, idsEligStart)
prepElig[idsEligStart] <- 1

vecStart <- rbinom(length(idsEligStart), 1, prep.start.prob)
idsStart <- idsEligStart[which(vecStart == 1)]
PrEPStarts <- length(idsStart)

# Set attributes for starters
if (length(idsStart) > 0) {
prepStat[idsStart] <- 1
prepStartTime[idsStart] <- at
prepLastRisk[idsStart] <- at
}

# PrEP adherence class
needPC <- which(is.na(prepClass[idsStart]))
prepClass[idsStart[needPC]] <- sample(x = 1:3, size = length(needPC),
replace = TRUE, prob = prep.adhr.dist)
## PrEP adherence class
needPC <- which(is.na(prepClass[idsStart]))
prepClass[idsStart[needPC]] <- sample(x = 1:3, size = length(needPC),
replace = TRUE, prob = prep.adhr.dist)

## Initiation (intervention; assuming effects rollling)
OptimInitStarts <- 0
PrEPStartsOptim <- 0
if (at >= prep.optim.start) {
# Starts (assumes intervention eligible one time rather than rolling)
idsEligStartInit <- intersect(idsNewIndic,
which(prepStat == 0 & prepElig == 1 & is.na(pOptimInitStatus)))
idsStartInit <- idsEligStartInit[which(rbinom(length(idsEligStartInit), 1, prep.optim.init.prob) == 1)]
OptimInitStarts <- length(idsStartInit)
pOptimInitStatus[idsStartInit] <- 1

# Stops
idsEligStopInit <- which(!is.na(pOptimInitStatus) & prepElig == 0)
pOptimInitStatus[idsEligStopInit] <- NA

# PrEP initiation on intervention
idsEligStartOptim <- intersect(idsIndic,
which(prepStat == 0 & pOptimInitStatus == 1))

vecStartOptim <- rbinom(length(idsEligStartOptim), 1, prep.start.prob.optim)
idsStartOptim <- idsEligStartOptim[which(vecStartOptim == 1)]
PrEPStartsOptim <- length(idsStartOptim)

# Set attributes for starters
if (length(idsStartOptim) > 0) {
prepStat[idsStartOptim] <- 1
prepStartTime[idsStartOptim] <- at
prepLastRisk[idsStartOptim] <- at
}

## PrEP adherence class
needPCoptim <- which(is.na(prepClass[idsStartOptim]))
prepClass[idsStartOptim[needPCoptim]] <- sample(x = 1:3, size = length(needPCoptim),
replace = TRUE, prob = prep.adhr.dist)
}


# Adherence & Retention Interventions -------------------------------------

OptimAdhrStarts <- 0
OptimRetnStarts <- 0
if (at >= prep.optim.start) {
## Adherence
# Starts (assumes intervention eligible one time rather than rolling)
idsEligStartAdhr <- which(prepStat == 1 & prepStartTime == at)
adhrSlotsOpen <- round(prep.optim.adhr.cap)
if (length(idsEligStartAdhr) > 0 & adhrSlotsOpen > 0) {
idsStartAdhrOptim <- ssample(idsEligStartAdhr,
min(length(idsEligStartAdhr), adhrSlotsOpen))
OptimAdhrStarts <- length(idsStartAdhrOptim)
pOptimAdhrStatus[idsStartAdhrOptim] <- 1
}

## PrEP adherence class (intervention)
if (OptimAdhrStarts > 0) {
prepClass[idsStartAdhrOptim] <- sample(x = 1:3,
size = OptimAdhrStarts,
replace = TRUE,
prob = prep.adhr.dist.optim)
}

# Stops
idsEligStopAdhr <- intersect(idsStp, which(!is.na(pOptimAdhrStatus)))
pOptimAdhrStatus[idsEligStopAdhr] <- NA


## Retention
# Starts (assumes intervention eligible one time rather than rolling)
idsEligStartRetn <- which(prepStat == 1 & prepStartTime == at)
retnCurrCov <- sum(pOptimRetnStatus == 1, na.rm = TRUE)
retnSlotsOpen <- round(prep.optim.retn.cap - retnCurrCov)
if (length(idsEligStartRetn) > 0 & retnSlotsOpen > 0) {
idsStartRetnOptim <- ssample(idsEligStartRetn,
min(length(idsEligStartRetn), retnSlotsOpen))
OptimRetnStarts <- length(idsStartRetnOptim)
pOptimRetnStatus[idsStartRetnOptim] <- 1
}

# Stops
idsEligStopRetn <- intersect(idsStp, which(!is.na(pOptimRetnStatus)))
pOptimRetnStatus[idsEligStopRetn] <- NA
}


Expand All @@ -163,6 +274,25 @@ prep_msm <- function(dat, at) {
dat$attr$prepLastRisk <- prepLastRisk
dat$attr$prepLastStiScreen <- prepLastStiScreen

dat$attr$pOptimInitStatus <- pOptimInitStatus
dat$attr$pOptimAdhrStatus <- pOptimAdhrStatus
dat$attr$pOptimRetnStatus <- pOptimRetnStatus

# Summary Stats
dat$epi$OptimInitStarts[at] <- OptimInitStarts
dat$epi$OptimInitPrev[at] <- sum(pOptimInitStatus == 1, na.rm = TRUE)
dat$epi$PrEPStarts[at] <- PrEPStarts
dat$epi$PrEPStartsOptim[at] <- PrEPStartsOptim
dat$epi$OptimAdhrStarts[at] <- OptimAdhrStarts
dat$epi$OptimAdhrPrev[at] <- sum(pOptimAdhrStatus == 1, na.rm = TRUE)
dat$epi$OptimRetnStarts[at] <- OptimRetnStarts
dat$epi$OptimRetnPrev[at] <- sum(pOptimRetnStatus == 1, na.rm = TRUE)

dat$epi$PrEPHighAdr[at] <- mean(dat$attr$prepClass == 3, na.rm = TRUE)
dat$epi$PrEPStopsInd[at] <- PrEPStopsInd
dat$epi$PrEPStopsRand[at] <- PrEPStopsRand
dat$epi$PrEPStopsRandOptim[at] <- PrEPStopsRandOptim

return(dat)
}

Expand Down
17 changes: 0 additions & 17 deletions R/mod.prevalence.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,23 +178,6 @@ prevalence_msm <- function(dat, at) {
dat$epi$cc.dx.delay.int.W[at] <- mean(diag.time[diag.time >= 3380 & race == 3] -
inf.time[diag.time >= 3380 & race == 3], na.rm = TRUE)

# same as above, but with medians
dat$epi$cc.dx.delay.med[at] <- median(diag.time[diag.time >= 2] - inf.time[diag.time >= 2], na.rm = TRUE)
dat$epi$cc.dx.delay.B.med[at] <- median(diag.time[diag.time >= 2 & race == 1] -
inf.time[diag.time >= 2 & race == 1], na.rm = TRUE)
dat$epi$cc.dx.delay.H.med[at] <- median(diag.time[diag.time >= 2 & race == 2] -
inf.time[diag.time >= 2 & race == 2], na.rm = TRUE)
dat$epi$cc.dx.delay.W.med[at] <- median(diag.time[diag.time >= 2 & race == 3] -
inf.time[diag.time >= 2 & race == 3], na.rm = TRUE)

dat$epi$cc.dx.delay.int.med[at] <- median(diag.time[diag.time >= 3380] - inf.time[diag.time >= 3380], na.rm = TRUE)
dat$epi$cc.dx.delay.int.B.med[at] <- median(diag.time[diag.time >= 3380 & race == 1] -
inf.time[diag.time >= 3380 & race == 1], na.rm = TRUE)
dat$epi$cc.dx.delay.int.H.med[at] <- median(diag.time[diag.time >= 3380 & race == 2] -
inf.time[diag.time >= 3380 & race == 2], na.rm = TRUE)
dat$epi$cc.dx.delay.int.W.med[at] <- median(diag.time[diag.time >= 3380 & race == 3] -
inf.time[diag.time >= 3380 & race == 3], na.rm = TRUE)

# dat$epi$cc.tx.any1y[at] <- sum((at - dat$attr$tx.period.last <= 52), na.rm = TRUE) /
# sum(dat$attr$diag.status == 1, na.rm = TRUE)
# dat$epi$cc.tx.any1y.B[at] <- sum((at - dat$attr$tx.period.last <= 52) & race == 1, na.rm = TRUE) /
Expand Down
26 changes: 13 additions & 13 deletions R/params.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -172,19 +172,19 @@
param_msm <- function(netstats,

# Clinical
hiv.test.rate = c(0.01325, 0.0125, 0.0124),
hiv.test.late.prob = c(0.25, 0.25, 0.25),
hiv.test.rate = c(0.00385, 0.00380, 0.00690),
hiv.test.late.prob = c(0, 0, 0),
test.window.int = 21/7,
tt.part.supp = c(0.20, 0.20, 0.20),
tt.full.supp = c(0.40, 0.40, 0.40),
tt.dur.supp = c(0.40, 0.40, 0.40),
tx.init.prob = c(0.092, 0.092, 0.127),
tx.halt.part.prob = c(0.0102, 0.0102, 0.0071),
tx.halt.full.rr = c(0.9, 0.9, 0.9),
tx.halt.dur.rr = c(0.5, 0.5, 0.5),
tx.reinit.part.prob = c(0.00066, 0.00066, 0.00291),
tx.reinit.full.rr = c(1.0, 1.0, 1.0),
tx.reinit.dur.rr = c(1.0, 1.0, 1.0),
tt.part.supp = c(0, 0, 0),
tt.full.supp = c(1, 1, 1),
tt.dur.supp = c(0, 0, 0),
tx.init.prob = c(0.1775, 0.190, 0.2521),
tx.halt.part.prob = c(0.0062, 0.0055, 0.0031),
tx.halt.full.rr = c(0.45, 0.45, 0.45),
tx.halt.dur.rr = c(0.45, 0.45, 0.45),
tx.reinit.part.prob = c(0.00255, 0.00255, 0.00255),
tx.reinit.full.rr = c(1, 1, 1),
tx.reinit.dur.rr = c(1, 1, 1),

# HIV natural history
max.time.off.tx.full.int = 52 * 15,
Expand All @@ -202,7 +202,7 @@ param_msm <- function(netstats,
vl.tx.down.slope = 0.25,
vl.tx.aids.down.slope = 0.25,
vl.tx.up.slope = 0.25,
aids.mr = 1/104,
aids.mr = 1/250,

# Demographic
a.rate = 0.00052,
Expand Down
Loading

0 comments on commit 962e109

Please sign in to comment.