## ----echo = FALSE, include = FALSE--------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----echo=FALSE, message=FALSE------------------------------------------------
library(EpiModel)

## ----setup, results = "hide"--------------------------------------------------
set.seed(10)

nw <- network_initialize(n = 200)
est <- netest(nw,
  formation = ~edges, target.stats = 60,
  coef.diss = dissolution_coefs(~offset(edges), 10, 0),  # duration = 10, closed population
  verbose = FALSE
)

param <- param.net(inf.prob = 0.9, rec.rate = 0.01, act.rate = 2)
init <- init.net(i.num = 10)
control <- control.net(type = "SIS", nsims = 1, nsteps = 250, verbose = FALSE)

## ----scenarios_setup----------------------------------------------------------
library(dplyr)

scenarios.df <- tribble(
  ~.scenario.id, ~.at, ~inf.prob, ~rec.rate,
  "base", 0, 0.9, 0.01,
  "initial_change", 0, 0.2, 0.01,
  "multiple_changes", 0, 0.1, 0.04,
  "multiple_changes", 100, 0.9, 0.01,
  "multiple_changes", 200, 0.1, 0.1
)

knitr::kable(scenarios.df)

## ----scenarios_list-----------------------------------------------------------
scenarios.list <- create_scenario_list(scenarios.df)
str(scenarios.list, max.level = 2)

## ----simulations--------------------------------------------------------------
# List to hold simulation results
d_list <- vector(mode = "list", length = length(scenarios.list))
names(d_list) <- names(scenarios.list)

for (scenario in scenarios.list) {
  print(scenario$id)
  sc.param <- use_scenario(param, scenario)
  sim <- netsim(est, sc.param, init, control)
  d_sim <- as.data.frame(sim)
  d_sim[["scenario"]] <- scenario$id
  d_list[[scenario$id]] <- d_sim
}

## ----plotting, fig.width = 8--------------------------------------------------
plot(d_list$base$time, d_list$base$i.num,
     type = "l", col = 1, lwd = 2, ylim = c(0, 250),
     xlab = "Time Step", ylab = "Number Infected")
lines(d_list$initial_change$time, d_list$initial_change$i.num,
      type = "l", col = 2, lwd = 2)
lines(d_list$multiple_changes$time, d_list$multiple_changes$i.num,
      type = "l", col = 3, lwd = 2)
abline(v = c(100, 200), lty = 2)
legend("topleft", legend = names(d_list),
       col = 1:3, lwd = 2, cex = 0.9, bty = "n")

## ----multi_param--------------------------------------------------------------
scenarios.df <- tribble(
  ~.scenario.id, ~.at, ~hiv.test.rate_1, ~hiv.test.rate_2, ~hiv.test.rate_3,
  "base", 0, 0.001, 0.001, 0.001,
  "initial_change", 0, 0.002, 0.001, 0.002,
  "multiple_changes", 0, 0.002, 0.001, 0.002,
  "multiple_changes", 100, 0.004, 0.001, 0.004,
  "multiple_changes", 200, 0.008, 0.001, 0.008
)
knitr::kable(scenarios.df)

## ----eval = FALSE-------------------------------------------------------------
# # CSV
# write.csv(scenarios.df, "scenarios.csv", row.names = FALSE)
# scenarios.df <- read.csv("scenarios.csv")
# 
# # RDS
# saveRDS(scenarios.df, "scenarios.rds")
# scenarios.df <- readRDS("scenarios.rds")

## ----param_table--------------------------------------------------------------
df_params <- tribble(
  ~param, ~value, ~type,
  "hiv.test.rate_1",  "0.003",        "numeric",
  "hiv.test.rate_2",  "0.102",        "numeric",
  "hiv.test.rate_3",  "0.492",        "numeric",
  "prep.require.lnt", "TRUE",         "logical",
  "group_1",          "first",        "character",
  "group_2",          "second",       "character"
)
knitr::kable(df_params)

## ----param_table2-------------------------------------------------------------
param <- param.net(data.frame.params = df_params,
                   other.param = c(5, 10), act.rate = 1)
param

## ----uncertainty_draws--------------------------------------------------------
set.seed(12)
n.draws <- 8

draws <- data.frame(
  .scenario.id = paste0("draw", seq_len(n.draws)),
  .at = 0,
  act.rate = sample(c(1, 2, 3), n.draws, replace = TRUE),
  inf.prob = rbeta(n.draws, 2, 6)
)
knitr::kable(draws, digits = 3)

## ----uncertainty_posterior----------------------------------------------------
posterior <- data.frame(act.rate = runif(500, 1, 3))
posterior$inf.prob <- 0.5 / posterior$act.rate * exp(rnorm(500, 0, 0.1))
cor(posterior$act.rate, posterior$inf.prob)

## ----uncertainty_joint--------------------------------------------------------
draws.joint <- data.frame(
  .scenario.id = paste0("draw", seq_len(n.draws)),
  .at = 0,
  posterior[sample(nrow(posterior), n.draws), ],
  row.names = NULL
)
knitr::kable(draws.joint, digits = 3)

## ----uncertainty_vector-------------------------------------------------------
draws.vec <- data.frame(
  .scenario.id = paste0("draw", 1:3),
  .at = 0,
  test.rate_1 = rnorm(3, 0.05, 0.01),
  test.rate_2 = rnorm(3, 0.15, 0.03)
)
knitr::kable(draws.vec, digits = 3)

param.vec <- param.net(inf.prob = 0.3, test.rate = c(0.05, 0.15))
sc.vec <- create_scenario_list(draws.vec)
use_scenario(param.vec, sc.vec$draw2)$test.rate

## ----uncertainty_run----------------------------------------------------------
param <- param.net(inf.prob = 0.3, act.rate = 2, rec.rate = 0.05)
control <- control.net(type = "SIS", nsims = 2, nsteps = 100, verbose = FALSE)

draws.list <- create_scenario_list(draws)

d_list <- lapply(draws.list, function(scenario) {
  sc.param <- use_scenario(param, scenario)
  sim <- netsim(est, sc.param, init, control)
  d_sim <- as.data.frame(sim)
  d_sim[["scenario"]] <- scenario$id
  d_sim
})
d_unc <- bind_rows(d_list)

## ----uncertainty_join---------------------------------------------------------
d_unc <- left_join(d_unc, select(draws, -.at), by = c("scenario" = ".scenario.id"))
d_unc$trans.prob <- 1 - (1 - d_unc$inf.prob)^d_unc$act.rate

d_final <- filter(d_unc, time == max(time))

## ----uncertainty_plot, fig.width = 8------------------------------------------
plot(d_final$trans.prob, d_final$i.num / d_final$num,
     pch = 19, col = d_final$act.rate, ylim = c(0, 1),
     xlab = "Transmission probability per partnership per time step",
     ylab = "Prevalence at time step 100")
legend("topleft", legend = paste("act.rate =", 1:3),
       col = 1:3, pch = 19, bty = "n")

## ----updater-example----------------------------------------------------------
list(
  at = 10,
  param = list(
    inf.prob = 0.3,
    act.rate = 0.5
  )
)

## ----updater-list-------------------------------------------------------------
list.of.updaters <- list(
  list(
    at = 100,
    param = list(
      inf.prob = 0.3,
      act.rate = 0.3
    )
  ),
  list(
    at = 125,
    param = list(
      inf.prob = 0.01
    )
  )
)

## ----updater-module-----------------------------------------------------------
param <- param.net(
 inf.prob = 0.1,
 act.rate = 0.1,
 .param.updater.list = list.of.updaters
)
init <- init.net(i.num = 10)
control <- control.net(
 type = "SI",
 nsims = 1,
 nsteps = 200,
 verbose = FALSE
)

## ----updater-module2----------------------------------------------------------
nw <- network_initialize(n = 100)
est <- netest(
  nw,
  formation = ~edges,
  target.stats = 50,
  coef.diss = dissolution_coefs(~offset(edges), 10, 0),  # duration = 10, closed population
  verbose = FALSE
)
mod <- netsim(est, param, init, control)

## -----------------------------------------------------------------------------
plot(mod, mean.smooth = FALSE)
abline(v = c(100, 125), lty = 2, col = "grey50")

## ----updater-verbose----------------------------------------------------------
list(
  at = 10,
  param = list(
    inf.prob = 0.3,
    act.rate = 0.5
  ),
  verbose = TRUE
)

## ----updater-function---------------------------------------------------------
list(
  at = 10,
  param = list(
    inf.prob = function(x) plogis(qlogis(x) + log(2)),
    act.rate = 0.5
  )
)

## ----control-updater-example--------------------------------------------------
list.of.updaters <- list(
  list(
    at = 100,
    control = list(
      resimulate.network = FALSE
    )
  ),
  list(
    at = 125,
    control = list(
      verbose = FALSE
    )
  )
)

## -----------------------------------------------------------------------------
control <- control.net(
 type = "SI",
 nsims = 1,
 nsteps = 200,
 verbose = TRUE,
 .control.updater.list = list.of.updaters
)