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

Scenario 1: Networks at equilibrium

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

# load files
files <- list.files("data/output",
  pattern = "(default)|(percent)|(global)", # nopatho",
  full.names = T
)

Potentially restrict which outputs to process.

# list parameter combinations and unique output files
param_files <- list.files(
  "data/parameters",
  pattern = "csv", full.names = T
)
params <- lapply(param_files, fread)
tags <- str_extract(param_files, pattern = "_[^_]+$") |>
  str_remove(pattern = c("_")) |>
  str_remove(pattern = ".csv")

# add scenario tags
params <- Map(params, tags, f = function(df, name) {
  df$scenario_tag <- name
  df
})

params <- rbindlist(params, use.names = TRUE, fill = TRUE)

# subset files for default scenario
params <- params[scenario_tag %in% c("default", "percent", "global"), ]
params <- params[(scenario_tag == "default" & costInfect == 0.25 &
  regen_time == 50) |
  (scenario_tag == "global" & costInfect == 0.25 &
    regen_time == 50) |
  (scenario_tag == "percent" & costInfect == 0.05 &
    regen_time == 50), ]

uids <- str_extract(files, "\\d{10}")
files <- files[uids %in% as.character(params$seed)]

Process network data

networks_folder <- "data/results/networks"
if (!dir.exists(networks_folder)) {
  message("Networks folder missing; creating")
  dir.create(networks_folder, recursive = TRUE)
}

for (i in files) {
  uid <- str_extract(i, "\\d{10}")

  output <- readRDS(i)[["output"]]

  nt <- pathomove::get_networks(output = output, assoc_threshold = 1)

  saveRDS(nt, file = glue("{networks_folder}/data_networks_{uid}.Rds"))
}

Quantifying transmission chains

# introduction generation
sgen <- 3000
generations <- c(seq(sgen, sgen + 100), seq(sgen + 500, sgen + 600)) # generation of introduction, + 200

if (!dir.exists("data/results/transmission_chains")) {
  dir.create("data/results/transmission_chains")
}

Save transmission chains

lapply(files, function(file) {
  output <- readRDS(file)

  gen_data <- pathomove::get_trait_data(output$output)

  # scenario tag
  scenario_tag <- output$scenario_tag

  gen_data <- gen_data[gen %in% generations, ]
  pathomove::get_social_strategy(gen_data)
  gen_data <- split(
    gen_data,
    by = "gen"
  )

  ## get transmission chains for visualisation
  transmission_chains <- lapply(
    gen_data[c("3000", "3500")], pathomove::get_transmission_chain
  )

  saveRDS(
    transmission_chains,
    file = sprintf(
      "data/results/transmission_chains/data_chains_%s_%i.Rds",
      output$scenario_tag, output$params$seed
    )
  )
})

# process secondary infection data and summarise
transmission_data <- lapply(files, function(file) {
  output <- readRDS(file)

  gen_data <- pathomove::get_trait_data(output$output)

  # scenario tag
  scenario_tag <- output$scenario_tag

  gen_data <- gen_data[gen %in% generations, ]
  pathomove::get_social_strategy(gen_data)
  gen_data <- split(
    gen_data,
    by = "gen"
  )

  # bind gen data again
  gen_data <- rbindlist(gen_data)

  # assign an id to all individuals
  gen_data[, id := seq(500), by = "gen"]

  # get all infected individuals, including those infected in the final timestep
  infected <- gen_data[t_infec > 0 | (!is.na(src_infect) & (src_infect != 0)), ]

  # get all infection source data
  infect_source <- gen_data[!is.na(src_infect) & (src_infect != 0),
    list(n_infected = .N),
    by = c("src_infect", "gen")
  ]

  # combine infection source and infected agent data
  infection_data <- merge(infected, infect_source,
    by.x = c("id", "gen"),
    by.y = c("src_infect", "gen"), all = TRUE
  )

  # fill zeros for NAs in n_infected
  infection_data[, n_infected := nafill(n_infected, fill = 0)]

  # count forward infections or secondary cases
  infection_data <- infection_data[, list(freq = .N),
    by = c("gen", "n_infected", "social_strat")
  ]

  # get simulation stage
  infection_data[, stage := fifelse(gen <= sgen + 500, "pre", "post")]

  # summarise mean secondary infections
  infection_data <- infection_data[, list(
    mean_freq = round(mean(freq))
  ), by = c("n_infected", "stage", "social_strat")]

  # prepare parameters
  params <- output$params[names(output$params) != "seed"]
  infection_data[, names(params) := params]
  infection_data[, scenario_tag := scenario_tag]

  # return infection data
  infection_data
})

# collect all transmission chain data
transmission_data <- rbindlist(transmission_data, fill = TRUE)

# save data - replicates are separate!
fwrite(
  transmission_data,
  file = "data/results/transmission_data.csv"
)

ggplot(transmission_data[scenario_tag == "default"]) +
  stat_summary(
    aes(n_infected, mean_freq, fill = social_strat),
    # col = "black",
    geom = "col", position = "dodge"
  ) +
  facet_grid(stage ~ scenario_tag)