Skip to content

Code specific to tergmLite within EpiModelHIV

Jump to bottom Edit New page
Samuel Jenness edited this page Jan 6, 2017 · 6 revisions

tergmLite changes within EpiModelHIV

Most of these changes include references to dat$el or dat$p, which are the new list-based formats for the network

Acts module

Creating edgelist object with references to dat$el

## Lines 36-60 (Lines 43, 51, 59)
# Parameters
    ai.scale <- dat$param$ai.scale
    if (type == "main") {
      base.ai.BB.rate <- dat$param$base.ai.main.BB.rate
      base.ai.BW.rate <- dat$param$base.ai.main.BW.rate
      base.ai.WW.rate <- dat$param$base.ai.main.WW.rate
      fixed <- FALSE
      ptype <- 1
      el <- dat$el[[1]]
    }
    if (type == "pers") {
      base.ai.BB.rate <- dat$param$base.ai.pers.BB.rate
      base.ai.BW.rate <- dat$param$base.ai.pers.BW.rate
      base.ai.WW.rate <- dat$param$base.ai.pers.WW.rate
      fixed <- FALSE
      ptype <- 2
      el <- dat$el[[2]]
    }
    if (type == "inst") {
      base.ai.BB.rate <- 1
      base.ai.BW.rate <- 1
      base.ai.WW.rate <- 1
      fixed <- ifelse(ai.scale != 1, FALSE, TRUE)
      ptype <- 3
      el <- dat$el[[3]]
    }

Births module

  1. Updating dat$el for MSM and heterosexual populations using add_vertices function within tergmLite
## Lines 55-60 (Line 58)
# Update Networks
  if (nBirths > 0) {
    for (i in 1:3) {
      dat$el[[i]] <- tergmLite::add_vertices(dat$el[[i]], nBirths)
    }
  }

## Lines 206-210 (Line 209)
# Update Population Structure
  if (nBirths > 0) {
    dat <- setBirthAttr_het(dat, at, nBirths)
    dat$el <- tergmLite::add_vertices(dat$el, nBirths)
  }
  1. stopifnot check for dat$el length
## Lines 212-214 (Line 212)
  if (unique(sapply(dat$attr, length)) != attributes(dat$el)$n) {
    stop("mismatch between el and attr length in births mod")
  }

Deaths module

  1. Updating dat$el for MSM and heterosexual populations using delete_vertices function within tergmLite
## Lines 54-62 (Lines 57, 60)
  if (length(dth.all) > 0) {
    dat$attr$active[dth.all] <- 0
    for (i in 1:3) {
      dat$el[[i]] <- tergmLite::delete_vertices(dat$el[[i]], dth.all)
    }
    dat$attr <- deleteAttr(dat$attr, dth.all)
    if (unique(sapply(dat$attr, length)) != attributes(dat$el[[1]])$n) {
      stop("mismatch between el and attr length in death mod")
    }

## Lines 161-164 (Line 163)
 ## 4. Update Population Structure ##
  inactive <- which(dat$attr$active == 0)
  dat$el <- tergmLite::delete_vertices(dat$el, inactive)
  dat$attr <- deleteAttr(dat$attr, inactive)
  1. stopifnot check for dat$el length
## Lines 166-168 (Line 166)
  if (unique(sapply(dat$attr, length)) != attributes(dat$el)$n) {
    stop("mismatch between el and attr length in death mod")
  }

Disclosure module

  1. Creating edgelist object with references to dat$el
## Lines 36-61 (Lines 44, 54, 60)
# Parameters and network
    if (type == "main") {
      disc.outset.B.prob <- dat$param$disc.outset.main.B.prob
      disc.at.diag.B.prob <- dat$param$disc.at.diag.main.B.prob
      disc.post.diag.B.prob <- dat$param$disc.post.diag.main.B.prob
      disc.outset.W.prob <- dat$param$disc.outset.main.W.prob
      disc.at.diag.W.prob <- dat$param$disc.at.diag.main.W.prob
      disc.post.diag.W.prob <- dat$param$disc.post.diag.main.W.prob
      el <- dat$el[[1]]
    }

    if (type == "pers") {
      disc.outset.B.prob <- dat$param$disc.outset.pers.B.prob
      disc.at.diag.B.prob <- dat$param$disc.at.diag.pers.B.prob
      disc.post.diag.B.prob <- dat$param$disc.post.diag.pers.B.prob
      disc.outset.W.prob <- dat$param$disc.outset.pers.W.prob
      disc.at.diag.W.prob <- dat$param$disc.at.diag.pers.W.prob
      disc.post.diag.W.prob <- dat$param$disc.post.diag.pers.W.prob
      el <- dat$el[[2]]
    }

    if (type == "inst") {
      disc.inst.B.prob <- dat$param$disc.inst.B.prob
      disc.inst.W.prob <- dat$param$disc.inst.W.prob
      el <- dat$el[[3]]
    }
  1. Creating master edgelist using dat$el
## Lines 136-144 (Line 138)
  if (at > 2) {
    discl.list <- dat$temp$discl.list
    master.el <- rbind(dat$el[[1]], dat$el[[2]], dat$el[[3]])
    m <- which(match(discl.list[, 1] * 1e7 + discl.list[, 2],
                     uid[master.el[, 1]] * 1e7 + uid[master.el[, 2]]) |
               match(discl.list[, 2] * 1e7 + discl.list[, 1],
                     uid[master.el[, 1]] * 1e7 + uid[master.el[, 2]]))
    dat$temp$discl.list <- discl.list[m, ]
  }

Initialize module

  1. Create dat$el and dat$p for MSM, use stergm_prep and ergm_prep functions within tergmLite
### Lines 45-60
  dat$el <- list()
  dat$p <- list()
  for (i in 1:2) {
    dat$el[[i]] <- as.edgelist(nw[[i]])
    attributes(dat$el[[i]])$vnames <- NULL
    p <- tergmLite::stergm_prep(nw[[i]], x[[i]]$formation, x[[i]]$coef.diss$dissolution,
                                x[[i]]$coef.form, x[[i]]$coef.diss$coef.adj, x[[i]]$constraints)
    p$model.form$formula <- NULL
    p$model.diss$formula <- NULL
    dat$p[[i]] <- p
  }
  dat$el[[3]] <- as.edgelist(nw[[3]])
  attributes(dat$el[[3]])$vnames <- NULL
  p <- tergmLite::ergm_prep(nw[[3]], x[[3]]$formation, x[[3]]$coef.form, x[[3]]$constraints)
  p$model.form$formula <- NULL
  dat$p[[3]] <- p
  1. Update dat$el and dat$p for MSM population
## Lines 742-743
  dat$el <- x$el[[s]]
  dat$p <- x$p[[s]]
  1. Create dat$el and dat$p for heterosexual population
## Lines 789-797
 nw <- simulate(x$fit, control = control.simulate.ergm(MCMC.burnin = 1e6))

  dat$el <- as.edgelist(nw)
  attributes(dat$el)$vnames <- NULL
  p <- tergmLite::stergm_prep(nw, x$formation, x$coef.diss$dissolution, x$coef.form,
                              x$coef.diss$coef.adj, x$constraints)
  p$model.form$formula <- NULL
  p$model.diss$formula <- NULL
  dat$p <- p
  1. Reinitialize heterosexual network
## Lines 875-893
reinit_het <- function(x, param, init, control, s) {
  dat <- list()
  dat$el <- x$el[[s]]
  dat$param <- param
  dat$param$modes <- 1
  dat$control <- control
  dat$nwparam <- x$nwparam
  dat$epi <- sapply(x$epi, function(var) var[s])
  names(dat$epi) <- names(x$epi)
  dat$attr <- x$attr[[s]]
  dat$stats <- list()
  dat$stats$nwstats <- x$stats$nwstats[[s]]
  dat$temp <- list()

  dat$param$modes <- 1
  class(dat) <- "dat"

  return(dat)
}

Risk history module

  1. Get degree distribution values from dat$el
## Lines 46-48
  main.deg <- get_degree(dat$el[[1]])
  casl.deg <- get_degree(dat$el[[2]])
  inst.deg <- get_degree(dat$el[[3]])

Network resimulation module

  1. Get degree, update terms in edgelist, and create new edges
## Lines 23-91
  ## Main network
  nwparam.m <- EpiModel::get_nwparam(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]],
                                             coef.form = nwparam.m$coef.form,
                                             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)
  } else {
    new.edges.m <- attributes(dat$el[[1]])$changes
    new.edges.m <- new.edges.m[new.edges.m[, "to"] == 1, 1:2, drop = FALSE]
  }
  dat$temp$new.edges <- matrix(dat$attr$uid[new.edges.m], ncol = 2)


  ## Casual network
  nwparam.p <- EpiModel::get_nwparam(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]],
                                             coef.form = nwparam.p$coef.form,
                                             coef.diss = nwparam.p$coef.diss$coef.adj,
                                             save.changes = TRUE)

  if (at == 2) {
    new.edges.p <- matrix(dat$el[[2]], ncol = 2)
  } else {
    new.edges.p <- attributes(dat$el[[2]])$changes
    new.edges.p <- new.edges.p[new.edges.p[, "to"] == 1, 1:2, drop = FALSE]
  }
  dat$temp$new.edges <- rbind(dat$temp$new.edges,
                              matrix(dat$attr$uid[new.edges.p], ncol = 2))


  ## One-off network
  nwparam.i <- EpiModel::get_nwparam(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]],
                                          coef = nwparam.i$coef.form)

  if (dat$control$save.nwstats == TRUE) {
    dat <- calc_resim_nwstats(dat, at)
  }

  return(dat)
}
  1. Calculate resimulated network statistics
## Lines 95-110 (Lines 98, 99, 101)
calc_resim_nwstats <- function(dat, at) {

  for (nw in 1:3) {
    n <- attr(dat$el[[nw]], "n")
    edges <- nrow(dat$el[[nw]])
    meandeg <- round(edges / n, 3)
    concurrent <- round(mean(get_degree(dat$el[[nw]]) > 1), 3)
    mat <- matrix(c(edges, meandeg, concurrent), ncol = 3, nrow = 1)
    if (at == 2) {
      dat$stats$nwstats[[nw]] <- mat
      colnames(dat$stats$nwstats[[nw]]) <- c("edges", "meand", "conc")
    }
    if (at > 2) {
      dat$stats$nwstats[[nw]] <- rbind(dat$stats$nwstats[[nw]], mat)
    }
  }
  1. Use edgelist to create attributes to update network parameters
## Lines 182-203
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")
    }
  1. Get degree distributions from dat$el
## Line 211
   dat$attr$deg.pers <- get_degree(dat$el[[2]])

## Line 289
   dat$attr$deg.pers <- get_degree(dat$el[[1]])

## Line 376 
   dat$attr$deg.pers <- get_degree(dat$el[[2]])
  1. Simulate edgelist for heterosexual population
## Lines 448-465
simnet_het <- function(dat, at) {

  # Update edges coefficients
  dat <- edges_correct_het(dat, at)

  # Update internal ergm data
  dat <- update_nwp_het(dat)

  # Pull network parameters
  nwparam <- get_nwparam(dat)

  # Simulate edgelist
  dat$el <- tergmLite::simulate_network(p = dat$p,
                                        el = dat$el,
                                        coef.form = nwparam$coef.form,
                                        coef.diss = nwparam$coef.diss$coef.adj)

  return(dat)
  1. Update network parameters for heterosexual population and update attributes
## Lines 468-475
update_nwp_het <- function(dat) {

  mf <- dat$p$model.form
  md <- dat$p$model.diss
  mhf <- dat$p$MHproposal.form
  mhd <- dat$p$MHproposal.diss

  n <- attributes(dat$el)$n

Transmission module

Create dat$el for transmission

## Lines 371-375
    if (is.null(dat$el)) {
      el <- get.dyads.active(dat$nw, at = at)
    } else {
      el <- dat$el
    }
Add a custom footer
Add a custom sidebar
Clone this wiki locally