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

Package index
Search the pratikunterwegs/patho-move-evol package
Vignettes
Functions
0
Source code
3
Man pages
0
Home
/
GitHub
/
pratikunterwegs/patho-move-evol
/

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

Plot evolutionary change across alternative implementations

library(data.table)
library(glue)

library(ggplot2)
library(patchwork)
library(colorspace)

Load data from default, percent, and global disperal scenarios

files <- list.files(
  "data/results/morph_data",
  pattern = "(default)|(percent)|(global)",
  full.names = TRUE
)

data_all <- lapply(files, fread)

# popsize and generation parameters
popsize <- 500
replicates <- 10
g_bin <- 100
g_increment <- 5 # data are for every 5th gen

# get time since pathogen
sgen <- 3000
genmax <- 5000

# regeneration rate
gen_time <- 100

# create columns for plotting
df_strat <- rbindlist(data_all)

# subset columns and rows to show 100th generation
df_strat[gen == max(gen), gen := genmax]
df_strat <- df_strat[gen %% 100 == 0 | gen == max(gen), ]
df_strat <- df_strat[, c(
  "gen", "social_strat", "N",
  "replicate", "scenario_tag", "costInfect",
  "regen_time", "dispersal"
)]

# calculate proportions
df_strat <- df_strat[, list(
  prop = sum(N) / (popsize * replicates)
), by = c(
  "gen", "social_strat",
  "scenario_tag", "costInfect",
  "regen_time", "dispersal"
)]

# express regeneration in generation time
df_strat[, regen_r := gen_time / regen_time]

df_strat <- split(
  df_strat,
  by = "scenario_tag"
)

Prepare strategy evolution plots for default, percent, and global scenarios

# prepare plots as a list
plots_evo <- lapply(
  df_strat,
  function(df) {
    p <-
      ggplot(df) +
      geom_col(
        aes(
          gen, prop,
          fill = social_strat
        ),
        width = 100,
        size = 1,
        position = "stack"
      ) +
      geom_vline(
        xintercept = sgen,
        lty = 2,
        linewidth = 0.3,
        col = "red"
      ) +
      scale_fill_discrete_sequential(
        palette = "Viridis",
        l2 = 80,
        rev = F,
        name = NULL,
        limits = c("agent avoiding", "agent tracking", "handler tracking"),
        labels = stringr::str_to_sentence,
        na.value = "lightgrey"
      ) +
      scale_x_continuous(
        breaks = c(1, sgen, genmax),
        name = "Generations",
        labels = scales::comma,
        sec.axis = dup_axis(
          breaks = sgen,
          labels = "Pathogen introduction",
          name = "Increasing productivity (social information less useful) \u279c",
        )
      ) +
      scale_y_continuous(
        labels = scales::percent,
        breaks = NULL,
        sec.axis = dup_axis(
          breaks = c(0, 0.5, 1),
          labels = scales::percent,
          name = "% Individuals"
        )
      ) +
      facet_grid(
        costInfect ~ regen_r,
        as.table = F,
        switch = c("y"),
        labeller = labeller(
          costInfect = function(x) {
            if (x >= 0.1) {
              sprintf("δE = %s", x)
            } else {
              scales::percent(as.numeric(x), prefix = "δE = ")
            }
          },
          regen_r = function(x) sprintf("R = %s times/gen", x)
        )
      ) +
      coord_cartesian(
        xlim = c(0, genmax),
        ylim = c(0, 1),
        expand = F
      ) +
      labs(
        y = glue::glue("Increasing infection cost \u279c")
      ) +
      theme_test(
        base_size = 8,
        base_family = "Arial"
      ) +
      theme(
        legend.position = "top",
        legend.key.height = unit(1, "mm"),
        strip.text = ggtext::element_markdown(),
        axis.text.x = element_text(hjust = 0.5, size = 6),
        axis.text.x.top = element_text(
          colour = "red"
        )
      )

    p
  }
)

# save fig 05
ggsave(
  plot = plots_evo[["default"]],
  filename = "figures/fig_evo_strategy_default.png",
  height = 120, width = 120, units = "mm"
)

# save supplementary figure for percent costs
ggsave(
  plot = plots_evo[["percent"]],
  filename = "supplement/figures/fig_evo_change_percent_cost.png",
  height = 120, width = 120, units = "mm"
)

# save supplementary figure for global dispersal
plots_evo$global <- plots_evo[["global"]] +
  labs(
    y = "Infection cost"
  )

Evolutionary outcomes in the threshold, handling time, and spatial scenarios

# read in alternative implementations
files_extra <- list.files(
  "data/results/morph_data",
  pattern = "(threshold)|(handling)|(spatial)",
  full.names = TRUE
)

# read in default scenario combinations as well
files_default <- list.files(
  "data/results/morph_data",
  pattern = "(default)",
  full.names = TRUE
)

# read in data
data_all <- lapply(files_extra, fread)
data_default <- lapply(files_default, fread)

# add data from default under the scenario_tag of "spatial"
df_default <- rbindlist(data_default)
df_default <- df_default[regen_time == 50 & costInfect == 0.25, ]
df_default[, scenario_tag := "spatial"]

# popsize and generation parameters
popsize <- 500
replicates <- 10
g_bin <- 100
g_increment <- 5 # data are for every 5th gen

# get time since pathogen
sgen <- 3000
genmax <- 5000

# regeneration rate
gen_time <- 100

# create columns for plotting
df_strat <- rbindlist(data_all)

# join with extra scenario data
df_strat <- rbindlist(list(df_strat, df_default), fill = TRUE)

# subset columns and rows to show 100th generation
df_strat[gen == max(gen), gen := genmax]
df_strat <- df_strat[gen %% 100 == 0 | gen == max(gen), ]
df_strat <- df_strat[, c(
  "gen", "social_strat", "N",
  "replicate", "scenario_tag", "costInfect",
  "regen_time", "nClusters", "clusterSpread", "handling_time"
)]

# calculate proportions
df_strat <- df_strat[, list(
  prop = sum(N) / (popsize * replicates)
), by = c(
  "gen", "social_strat",
  "scenario_tag", "costInfect",
  "regen_time", "nClusters", "handling_time", "clusterSpread"
)]

# edit number of food patches for uniform distribution case
df_strat$nClusters <- as.character(df_strat$nClusters)
df_strat[
  scenario_tag == "spatial" & nClusters == 10 & clusterSpread == 10,
  nClusters := "uniform"
]

# express regeneration in generation time
df_strat[, regen_r := gen_time / regen_time]

# factor levels for nCLusters
df_strat$nClusters <- factor(
  df_strat$nClusters,
  levels = c("uniform", "60", "10")
)

Prepare figures for minor scenarios

df_strat <- split(df_strat, by = "scenario_tag")
df_strat <- df_strat[c("threshold", "handling", "spatial")]

plot_list <- lapply(df_strat, function(df) {
  ggplot(df) +
    geom_col(
      aes(
        gen, prop,
        fill = social_strat
      ),
      width = 100,
      size = 1,
      position = "stack"
    ) +
    geom_vline(
      xintercept = sgen,
      lty = 2,
      linewidth = 0.3,
      col = "red"
    ) +
    scale_fill_discrete_sequential(
      palette = "Viridis",
      l2 = 80,
      rev = F,
      name = NULL,
      limits = c("agent avoiding", "agent tracking", "handler tracking"),
      labels = stringr::str_to_sentence,
      na.value = "lightgrey"
    ) +
    scale_x_continuous(
      breaks = c(1, sgen, genmax),
      labels = scales::comma,
      sec.axis = dup_axis(
        breaks = sgen,
        labels = "Pathogen introduction",
        name = NULL
      )
    ) +
    scale_y_continuous(
      labels = scales::percent,
      breaks = NULL,
      sec.axis = dup_axis(
        breaks = c(0, 0.5, 1),
        labels = scales::percent,
        name = "% Individuals"
      )
    ) +
    coord_cartesian(
      xlim = c(0, genmax),
      ylim = c(0, 1),
      expand = F
    )
})

# for scenario with reproduction threshold
plot_list$threshold <-
  plot_list$threshold +
  facet_grid(
    cols = vars(costInfect),
    rows = vars(regen_r),
    as.table = F,
    switch = c("y"),
    labeller = labeller(
      costInfect = function(x) {
        if (x >= 0.1) {
          sprintf("δE = %s", x)
        } else {
          scales::percent(as.numeric(x), prefix = "δE = ")
        }
      },
      regen_r = function(x) sprintf("R = %s times/gen", x)
    )
  ) +
  labs(
    x = "Generations",
    y = "Landscape productivity"
  ) +
  guides(
    x.sec = guide_axis(
      title = "Increasing disease cost \u279c"
    )
  )

# for scenario with different handling times
plot_list$handling <-
  plot_list$handling +
  facet_grid(
    cols = vars(handling_time),
    rows = vars(regen_r, costInfect),
    as.table = F,
    switch = c("y"),
    labeller = labeller(
      handling_time = function(x) {
        sprintf("T<sub>H</sub> = %s", x)
      },
      costInfect = function(x) {
        if (x >= 0.1) {
          sprintf("δE = %s", x)
        } else {
          scales::percent(as.numeric(x), prefix = "δE = ")
        }
      },
      regen_r = function(x) sprintf("R = %s times/gen", x),
      .multi_line = FALSE
    )
  ) +
  labs(
    x = "Generations",
    y = "Land. product., Infection cost"
  ) +
  guides(
    x.sec = guide_axis(
      title = "Increasing handling time (soc. info. available longer) \u279c"
    )
  )

# plot with different spatial configurations
plot_list$spatial <-
  plot_list$spatial +
  facet_grid(
    cols = vars(nClusters),
    rows = vars(regen_r, costInfect),
    as.table = F,
    switch = c("y"),
    labeller = labeller(
      handling_time = function(x) {
        sprintf("T<sub>H</sub> = %s", x)
      },
      costInfect = function(x) {
        if (x >= 0.1) {
          sprintf("δE = %s", x)
        } else {
          scales::percent(as.numeric(x), prefix = "δE = ")
        }
      },
      regen_r = function(x) sprintf("R = %s times/gen", x),
      nClusters = function(x) {
        ifelse(x == "uniform",
          "Uniform distribution",
          sprintf("%s food patches", x)
        )
      },
      .multi_line = FALSE
    )
  ) +
  labs(
    x = "Generations",
    y = "Land. product., Infection cost"
  ) +
  guides(
    x.sec = guide_axis(
      title = "Increasing landscape patchiness \u279c"
    )
  )

# add formatting
plot_list <- lapply(plot_list, function(pl) {
  pl +
    theme_test(
      base_size = 8,
      base_family = "Arial"
    )
})

Plot threshold and handling.

# arrange alternative scenario plots
plot_alternatives <- wrap_plots(
  plots_evo$global,
  plot_list$threshold,
  plot_list$handling,
  guides = "collect",
  design = "A\nB\nC"
) &
  plot_annotation(tag_levels = c("A")) &
  theme(
    plot.tag = element_text(
      face = "bold", size = 12
    ),
    legend.position = "top",
    legend.key.height = unit(1, "mm"),
    strip.text = ggtext::element_markdown(),
    panel.border = element_rect(fill = NA),
    axis.text.x.top = element_text(
      colour = "red"
    )
  )

ggsave(
  plot_alternatives,
  file = "supplement/figures/fig_eco_evo_alternatives.png",
  height = 150, width = 120, units = "mm"
)

Plot evolutionary outcomes for spatial structure

Prepare landscape examples

# generate landscapes
landscapes <- Map(
  c(10, 60, 10),
  c(10, 1, 1),
  c("uniform", "60", "10"),
  f = function(clusters, spread, name) {
    df <- pathomove::get_test_landscape(
      nItems = 1800,
      landsize = 60,
      nClusters = clusters, clusterSpread = spread,
      regen_time = 50
    )
    df$name <- name
    df
  }
)

landscapes <- rbindlist(landscapes)

# save landscapes
fwrite(
  landscapes,
  # file = "data/results/test_landscapes.csv" # commented for safety
)

landscapes <- fread("data/results/test_landscapes.csv")
landscapes$name <- factor(
  landscapes$name,
  levels = c("uniform", "60", "10")
)

Prepare landscape plots

plot_landscape <-
  ggplot(landscapes, aes(x, y)) +
  geom_point(
    size = 0.1,
    aes(col = tAvail),
    show.legend = FALSE
  ) +
  scale_colour_continuous_sequential(
    palette = "Greens 2",
    rev = FALSE
  ) +
  facet_grid(
    cols = vars(name),
    labeller = labeller(
      name = function(x) {
        ifelse(x == "uniform",
          "Uniform distribution",
          sprintf("%s food patches", x)
        )
      }
    )
  ) +
  coord_cartesian(
    expand = FALSE,
    xlim = c(0, 60),
    ylim = c(0, 60)
  ) +
  theme_test(base_family = "Arial", base_size = 8) +
  theme(
    strip.background = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank()
  )

Prepare outcomes and landscape plot

# arrange alternative scenario plots
plot_spatial <- wrap_plots(
  plot_list$spatial,
  plot_landscape,
  guides = "collect",
  design = "A\nB"
) &
  plot_annotation(tag_levels = c("A")) &
  theme(
    plot.tag = element_text(
      face = "bold", size = 12
    ),
    legend.position = "top",
    legend.key.height = unit(1, "mm"),
    strip.text = ggtext::element_markdown(),
    panel.border = element_rect(fill = NA),
    axis.text.x.top = element_text(
      colour = "red"
    )
  )

ggsave(
  plot_spatial,
  file = "supplement/figures/fig_eco_evo_spatial.png",
  height = 100, width = 120, units = "mm"
)

Compare movement, intake, and associations in main scenarios

# list files, read, and filter for default model options
files <- list.files(
  "data/results/gen_data",
  pattern = "(default)|(percent)|(global)",
  full.names = TRUE
)

data_all <- lapply(files, fread)
data_all <- rbindlist(data_all)

data_all <- data_all[scenario == 1, ]
data <- copy(data_all)

# get time since pathogen
sgen <- 3000
genmax <- 5000
gen_time <- 100

# plot regen as rate per generation timesteps
data[, regen_r := gen_time / regen_time]

# select a range of generations before and after
data <- data[between(gen, sgen - 500, sgen) |
  between(gen, sgen + 500, sgen + 1000)]

# label pre and post pathogen generations
data[, period := fifelse(
  gen <= sgen, "pre", "post"
)]

data[, period := factor(
  period,
  levels = c("pre", "post")
)]

# split by simulation type
data <- split(
  data,
  by = "scenario_tag"
)

Plot ecological outcomes in default and minor scenarios

plots_eco <- Map(
  data, names(data),
  f = function(df, nm) {
    df <- melt(
      df[, c(
        "intake.mean", "moved.mean", "assoc.mean",
        "period", "costInfect", "regen_r", "replicate"
      )],
      id.vars = c("period", "costInfect", "regen_r", "replicate")
    )

    df <- split(
      df,
      by = "variable"
    )

    df <- df[sprintf("%s.mean", c("moved", "intake", "assoc"))]

    col_pal <- colorspace::diverging_hcl(
      3,
      palette = "Tofino", l = 50, c = 80
    )[c(1, 3, 2)]

    names <- c("Mean distance moved", "Per-capita intake", "Mean associations")

    p <- Map(
      df, col_pal, names,
      f = function(le, col, n) {
        pl <- ggplot(le) +
          stat_summary(
            fun = median,
            fun.min = function(x) median(x) - sd(x),
            fun.max = function(x) median(x) + sd(x),
            aes(
              period, value
            ),
            shape = 21,
            fill = col
          ) +
          facet_grid(
            costInfect ~ regen_r,
            as.table = F,
            switch = c("both"),
            labeller = labeller(
              cost = function(x) {
                if (x >= 0.1) {
                  sprintf("δE = %s", x)
                } else {
                  scales::percent(as.numeric(x), prefix = "δE = ")
                }
              },
              regen_r = function(x) sprintf("R = %s times/gen", x)
            )
          ) +
          scale_x_discrete(
            labels = c(
              "Pre-.",
              "Post-."
            ),
            name = "Increasing productivity (social information less useful) \u279c"
          ) +
          scale_y_continuous(
            name = n
          ) +
          coord_cartesian(
            expand = T
          ) +
          theme_test(
            base_size = 8,
            base_family = "Arial"
          ) +
          theme(
            legend.position = "top",
            legend.key.height = unit(1, "mm"),
            strip.placement = "outside",
            strip.text = ggtext::element_markdown(),
            axis.text.x = element_text(hjust = 0.5, size = 6)
          )

        if (n == "Mean associations") {
          pl <- pl +
            scale_y_continuous(
              labels = scales::comma,
              name = "Mean associations"
            ) +
            theme(
              axis.text.y = element_text(
                angle = 90, hjust = 0.5
              )
            )
        }
        pl
      }
    )

    # wrap plots
    p <- wrap_plots(
      p
    ) &
      plot_annotation(
        tag_levels = "A"
      )
    # save data
    ggsave(
      p,
      filename = glue("supplement/figures/fig_eco_compare_{nm}.png"),
      height = 90,
      width = 240,
      units = "mm"
    )
  }
)


pratikunterwegs/patho-move-evol documentation built on April 15, 2023, 5:54 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close