In pratikunterwegs/patho-move-evol: Analysis code for the _pathomove_ model.

Run SIR models

Load libraries.

library(data.table)
library(igraph)
library(tidygraph)
library(glue)

Prepare parameters

# disease parameters
beta <- c(5)
gamma <- c(1)
threshold <- c(1, 10)

# make combinations
params_sir <- CJ(beta, gamma, threshold)

beta <- params_sir$beta
gamma <- params_sir$gamma
threshold <- params_sir$threshold

Load networks from default scenario

Read parameter files to subset for default parameter combination, based on the "ofi" column.

data_files <- list.files(
  "data/results/networks",
  full.names = TRUE,
  pattern = "Rds"
)

Run SIR models

sc_repl <- seq(10) # the number of replicates

# the generation of pathogen introduction
gen_patho_intro <- as.character(3000)
gen_patho_adapt <- as.character(3500)

data_sir_models <- Map(
  data_files, sc_repl,
  f = function(file, sc_repl) {

    # load networks
    ntwk <- readRDS(file)

    # pre disease
    d_pre <- Map(
      beta, gamma, threshold,
      f = function(b, g, thr) {
        nt <- ntwk[[gen_patho_intro]] # hardcoded for now

        params_ <- as.list(as_tibble(nt)[1, ])

        # filter for threshold
        nt <- nt %>%
          activate(edges) %>%
          filter(weight > thr) %>%
          activate(nodes)

        d_pre_ <- igraph::sir(
          graph = nt, beta = b, gamma = g, no.sim = 25
        ) |>
          pathomove::handle_sir_data(digits = 1)
        d_pre_$type <- "pre"

        # add parameters
        d_pre_[, names(params_) := params_]
        d_pre_[, c("beta", "gamma", "threshold") := list(b, g, thr)]

        d_pre_
      }
    ) |>
      rbindlist()

    # post disease
    d_post <- Map(
      beta, gamma, threshold,
      f = function(b, g, thr) {
        nt <- ntwk[[gen_patho_adapt]] # hardcoded for now

        params_ <- as.list(as_tibble(nt)[1, ])

        # filter for threshold
        nt <- nt %>%
          activate(edges) %>%
          filter(weight > thr) %>%
          activate(nodes)

        d_post_ <- igraph::sir(
          graph = nt, beta = b, gamma = g, no.sim = 25
        ) |>
          pathomove::handle_sir_data(digits = 1)
        d_post_$type <- "post"

        # add parameters
        d_post_[, names(params_) := params_]
        d_post_[, c("beta", "gamma", "threshold") := list(b, g, thr)]

        d_post_
      }
    ) |>
      rbindlist()

    data <- rbindlist(
      list(d_pre, d_post)
    )

    # add simulation replicate
    data$sc_repl <- sc_repl
    data
  }
)

data_sir_models <- rbindlist(data_sir_models)

# save data
fwrite(
  data_sir_models,
  file = "data/results/data_sir_models.csv"
)

# sanity check
data_sir_models[class != "NS" & time < 5] %>%
  ggplot() +
  stat_summary(
    aes(time, mean, col = type)
  ) +
  facet_grid(
    threshold ~ class,
    labeller = label_both
  ) +
  # scale_x_sqrt()+
  coord_cartesian(
    xlim = c(0, 5)
  )