_results/fig3.R

In lucasnell/evo_alt_states: Simulations for alternative stable states and coexistence in evolution

# start ----

# load packages and get helper functions
source("_results/_preamble.R")

# whether to re-do simulations (use rds files otherwise)
.REDO_SIMS <- FALSE
# whether to re-save plots
.RESAVE_PLOTS <- TRUE



# =============================================================================*
# =============================================================================*

# *--------------------------* ----
# Fig 3: Stochasticity ----

# =============================================================================*
# =============================================================================*


#'
#' .sigma_V should be wrapped into a list if you want it to vary by axis.
#'
one_stoch_sim <- function(.eta,
                          .sigma_V,
                          .d,
                          ...) {

    # .eta = 1e-2; .sigma_V = 0.05

    stopifnot(length(.d) == 2)
    stopifnot(sign(.d) == c(-1,1))

    # Used for all calls to `quant_gen` below
    args <- list(eta = .eta,
                 d = .d,
                 q = 2,
                 n = 4,
                 # These values were chosen because their distance from
                 # the origin is 1:
                 V0 = cbind(c(0.1961161, 0.9805805), c(0.6401844, 0.7682213),
                            c(0.9805805, 0.1961161), c(0.7682213, 0.6401844)),
                 sigma_V0 = 0,
                 final_t = 200e3L,
                 spp_gap_t = 0L,
                 save_every = 100L,
                 sigma_V = unlist(.sigma_V),
                 n_reps = 24,
                 n_threads = .N_THREADS,
                 show_progress = FALSE)

    others <- list(...)
    for (n in names(others)) args[[n]] <- others[[n]]
    args$n_threads <- min(args$n_reps, args$n_threads)

    nv <- args %>%
        do.call(what = quant_gen) %>%
        .$nv %>%
        pivot() %>%
        mutate(d1 = .d[1], d2 = .d[2]) %>%
        mutate(eta = .eta, sigma_V = .sigma_V) %>%
        select(eta, sigma_V, d1, d2, everything())

    return(nv)

}




# Values of sigma_V to use for each suite of traits:
sigma_v_vals <- c(0, 0.11, 0.2)
# To use for levels of vector combining sigmas in both suites of traits:
sigma_v_levels <- do.call(c, map(rev(sigma_v_vals),
                                 ~ paste(sigma_v_vals, .x, sep = "__")))
#
# z <- crossing(.eta = 1 * 1e-2,
#          .sigma_V1 = 0.05,
#          .sigma_V2 = 0.1,
#          .d = list(c(-1e-2, 1e-2))) %>%
#     mutate(.sigma_V = map2(.sigma_V1, .sigma_V2, ~ list(c(.x, .y)))) %>%
#     select(-.sigma_V1, -.sigma_V2) %>%
#     pmap_dfr(one_stoch_sim)
#
# z %>%
#     filter(time == max(time)) %>%
#     mutate(grp = group_spp(V1, V2, .prec = 0.1)) %>%
#     group_by(grp) %>%
#     summarize(n = n(),
#               V1 = mean(V1), V2 = mean(V2))





# Takes ~ 70 sec
set.seed(1909344759)
all_stoch_sims <- crossing(.eta = -1:1 * 1e-2,
                           .sigma_V1 = sigma_v_vals,
                           .sigma_V2 = sigma_v_vals,
                           .d = list(c(-1e-2, 1e-2))) %>%
    mutate(.sigma_V = map2(.sigma_V1, .sigma_V2, ~ list(c(.x, .y)))) %>%
    select(-.sigma_V1, -.sigma_V2) %>%
    pmap_dfr(one_stoch_sim) %>%
    mutate(sigma_v1 = map_dbl(sigma_V, ~ .x[1]),
           sigma_v2 = map_dbl(sigma_V, ~ .x[2]),
           sigma_V = paste(sigma_v1, sigma_v2, sep = "__")) %>%
    mutate(sigma_V = factor(sigma_V, levels = sigma_v_levels),
           eta = factor(sign(eta), levels = -1:1,
                        labels = c("Sub-additive", "Additive", "Super-additive"))) %>%
    select(-starts_with("Vp"), -d1, -d2)


# # saveRDS(all_stoch_sims, "./_results/results-rds/all_stoch_sims.rds")
#
# all_stoch_sims <- readRDS("./_results/results-rds/all_stoch_sims.rds")




#'
#' An alternative set of simulations for additive or strongly non-additive,
#' for the supplemental info.
#'

# Takes ~50 sec
set.seed(567843943)
all_stoch_sims_alt <- crossing(.eta = c(-1,1) * 0.6,
                                    .sigma_V1 = sigma_v_vals,
                                    .sigma_V2 = sigma_v_vals,
                                    .d = list(c(-1e-2, 1e-2))) %>%
    mutate(.sigma_V = map2(.sigma_V1, .sigma_V2, ~ list(c(.x, .y)))) %>%
    select(-.sigma_V1, -.sigma_V2) %>%
    pmap_dfr(one_stoch_sim) %>%
    mutate(sigma_v1 = map_dbl(sigma_V, ~ .x[1]),
           sigma_v2 = map_dbl(sigma_V, ~ .x[2]),
           sigma_V = paste(sigma_v1, sigma_v2, sep = "__")) %>%
    mutate(sigma_V = factor(sigma_V, levels = sigma_v_levels),
           eta = factor(sign(eta), levels = c(-1, 1),
                        labels = c("Sub-additive", "Super-additive"))) %>%
    select(-starts_with("Vp"), -d1, -d2)

# # saveRDS(all_stoch_sims_alt, "./_results/results-rds/all_stoch_sims_alt.rds")
#
# all_stoch_sims_alt <- readRDS("./_results/results-rds/all_stoch_sims_alt.rds")



# Shows that reps don't differ in their outcomes, so I'm using rep #1
# as a representative example from now on.
map(list(all_stoch_sims, all_stoch_sims_alt),
    function(.dd) {
        .dd %>%
            filter(time == max(time)) %>%
            group_by(eta, sigma_V) %>%
            # Group species by their ending investments:
            mutate(grp = group_spp(V1, V2, .prec = 0.1)) %>%
            group_by(eta, sigma_V, rep) %>%
            summarize(comm_conf = paste(sort(grp), collapse = "_"),
                      .groups = "drop") %>%
            group_by(eta, sigma_V) %>%
            # If comms differed by rep, then we'd see " xx " in the `comm_conf` col:
            summarize(comm_conf = paste(unique(comm_conf), collapse = " xx "),
                      .groups = "drop") %>%
            filter(grepl(" xx ", comm_conf))
    })

all_stoch_sims <- all_stoch_sims %>%
    filter(rep == 1) %>%
    select(-rep)
all_stoch_sims_alt <- all_stoch_sims_alt %>%
    filter(rep == 1) %>%
    select(-rep)





#'
#' Provides `X` and `Y` (or `Ymin` and `Ymax`) for plotting PDF of normal
#' distribution with a given mean and SD.
#' If `dens_on_x = FALSE` (the default), then it provides x values in the
#' `X` column and densities in the `Y` column, for use inside `geom_area`.
#' If `dens_on_x = TRUE`, then densities are on the x axis (so are in the `X`
#' column) and x values are in columns `Ymin` and `Ymax`;
#' this should be used in `geom_ribbon`.
#'
norm_pdf_data <- function(mean, sd, dens_on_x = FALSE, d_min = 1e-2) {
    if (sd <= 0) return(NULL)
    idn <- function(d,m) mean + m * sqrt(-2 * sd^2 *
                                             log(sd * sqrt(2 * pi) * d))
    if (dens_on_x) {
        d_max <- dnorm(mean, mean, sd)
        dens <- seq(d_min, d_max, length.out = 101)
        half1 <- sapply(dens, idn, m = -1)
        half2 <- sapply(dens, idn, m = 1)
        tibble(X = dens, Ymin = half1, Ymax = half2)
    } else {
        tibble(X = seq(idn(d_min, -1), idn(d_min, 1), length.out = 101)) %>%
            mutate(Y = map_dbl(X, ~ dnorm(.x, mean, sd)))
    }
}



#'
#' Plots how species evolve through axis space, for a given set of simulations
#' and `sigma_V` values.
#' Note that `.s` should be a single numeric or character.
#' In the latter case, it assumes that `.sims` combines multiple values
#' into a single factor with the separator `"__"`.
#' `.e` should be `"Super-additive"` or `"Sub-additive"` (case doesn't matter)
#' `.dens_mult` adjusts the height of the PDF curves
#'
stoch_comm_V_plotter <- function(.s, .e, .sims,
                                 .dens_mult = 0.075,
                                 theme_args = list()) {

    # .s = levels(all_stoch_sims$sigma_V)[5]; .e = "Super-additive"
    # .sims = all_stoch_sims; .dens_mult = 0.075; theme_args = list()
    # rm(.s, .e, .sims, .dens_mult, theme_args, .s_fct, .equil_comms,
    #    x_pdf, y_pdf, p)

    stopifnot(length(.s) == 1 && (is.numeric(.s) || is.character(.s)))
    if (is.character(.s)) {
        .s_fct <- .s
        .s <- as.numeric(str_split(.s, "__")[[1]])
    } else {
        .s_fct <- .s
        .s <- rep(.s, 2)
    }

    .e <- match.arg(str_to_sentence(.e), c("Super-additive", "Additive",
                                           "Sub-additive"))

    .sims_es <- .sims %>%
        filter(eta == .e, sigma_V == .s_fct) %>%
        select(-eta, -starts_with("sigma"), -N)

    .equil_comms <- .sims_es %>%
        filter(time == max(time)) %>%
        mutate(grp = group_spp(V1, V2, .prec = 1e-1)) %>%
        group_by(grp) %>%
        summarize(V1 = mean(V1), V2 = mean(V2), n_spp = n(),
                  .groups = "drop") %>%
        select(-grp) %>%
        arrange(V1, V2)

    # Filter out time points after equilibrium reached.
    min_asd <- function(.v1, .v2) {
        # This fxn calculates minimum average square differences
        .asd <- (.v1 - .equil_comms$V1)^2 + (.v2 - .equil_comms$V2)^2 / 2
        return(min(.asd))
    }
    for (t in sort(unique(.sims_es$time))) {
        .sims_es_t <- filter(.sims_es, time == t)
        masd <- map2_dbl(.sims_es_t$V1, .sims_es_t$V2, min_asd)
        if (max(masd) < 1e-4) break
    }
    if (t < max(.sims_es$time)) {
        .sims_es <- .sims_es %>% filter(time <= t)
    }

    if (.s[1] > 0) {
        x_pdf <- geom_area(data = norm_pdf_data(0.75, .s[1]),
                           aes(X, Y * .dens_mult), color = NA, fill = "gray80")
    } else x_pdf <- geom_blank()
    if (.s[2] > 0) {
        y_pdf <- geom_ribbon(data = norm_pdf_data(0.75, .s[2], TRUE),
                             aes(X * .dens_mult, ymin = Ymin, ymax = Ymax),
                             inherit.aes = FALSE, color = NA, fill = "gray80")
    } else y_pdf <- geom_blank()


    p <- .sims_es %>%
        ggplot(aes(V1, V2, color = spp)) +
        geom_abline(slope = 1, intercept = 0,
                    linetype = 2, color = "gray70") +
        x_pdf +
        y_pdf +
        geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
        geom_vline(xintercept = 0, color = "gray70", size = 0.5) +
        geom_path() +
        geom_point(data = filter(.sims_es, time == 0), shape = 1, size = 1) +
        geom_point(data = .equil_comms, shape = 19, color = "black", size = 3) +
        scale_color_viridis_d(option = "viridis", end = 0.9, begin = 0.4,
                              guide = "none") +
        scale_y_continuous("Partitioning investment", breaks = 0:1) +
        scale_x_continuous("Conflict investment", breaks = 0:1) +
        coord_equal(xlim = c(-0.1, 1.6), ylim = c(-0.1, 1.6)) +
        theme(panel.background = element_blank(),
              plot.margin = margin(0,0,0,0)) +
        do.call(theme, theme_args) +
        NULL

    return(p)

}



# ----------------------------------------*
# Main-text plot
# ----------------------------------------*

all_sub_p <- levels(all_stoch_sims$sigma_V) %>%
    map(~ stoch_comm_V_plotter(.x, "Sub-additive", all_stoch_sims)) %>%
    imap(adjust_grid_axes, rows = 3, cols = 3, no_titles = "xy") %>%
    wrap_plots(ncol = 3, byrow = TRUE)

all_add_p <- levels(all_stoch_sims$sigma_V) %>%
    map(~ stoch_comm_V_plotter(.x, "Additive", all_stoch_sims)) %>%
    imap(adjust_grid_axes, rows = 3, cols = 3, no_titles = "xy") %>%
    wrap_plots(ncol = 3, byrow = TRUE)

all_super_p <- levels(all_stoch_sims$sigma_V) %>%
    map(~ stoch_comm_V_plotter(.x, "Super-additive", all_stoch_sims)) %>%
    imap(adjust_grid_axes, rows = 3, cols = 3, no_titles = "xy") %>%
    wrap_plots(ncol = 3, byrow = TRUE)

all_stoch_layout <- "#246
                     1357
                     #888"

all_stoch_p_main <- wrap_elements(textGrob("Partitioning investment", rot = 90),
                                  clip = FALSE) +
    wrap_elements(textGrob("weak sub-additive", y = 1, vjust = 1), clip = FALSE) +
    all_sub_p +
    wrap_elements(textGrob("additive", y = 1, vjust = 1), clip = FALSE) +
    all_add_p +
    wrap_elements(textGrob("weak super-additive", y = 1, vjust = 1), clip = FALSE) +
    all_super_p +
    wrap_elements(textGrob("Conflict investment"), clip = FALSE) +
    plot_layout(heights = c(0.02, 1, 0.005),
                widths = c(0.01, 1, 1, 1),
                design = all_stoch_layout)

all_stoch_p <- function() {
    grid.newpage()
    grid.draw(all_stoch_p_main)
    grid.draw(textGrob("A", x = 0.05, y = 0.99, just = c("left", "top")))
    grid.draw(textGrob("B", x = 0.37, y = 0.99, just = c("left", "top")))
    grid.draw(textGrob("C", x = 0.69, y = 0.99, just = c("left", "top")))
    invisible(NULL)
}

# all_stoch_p()

# if (.RESAVE_PLOTS) save_plot(all_stoch_p, 9, 3.5, .prefix = "3-", .png = TRUE)



# ----------------------------------------*
# Supplement plot
# ----------------------------------------*


all_sub_alt_p <- levels(all_stoch_sims_alt$sigma_V) %>%
    map(~ stoch_comm_V_plotter(.x, "Sub-additive", all_stoch_sims_alt)) %>%
    imap(adjust_grid_axes, rows = 3, cols = 3, no_titles = "xy") %>%
    wrap_plots(ncol = 3, byrow = TRUE)

all_super_alt_p <- levels(all_stoch_sims_alt$sigma_V) %>%
    map(~ stoch_comm_V_plotter(.x, "Super-additive", all_stoch_sims_alt)) %>%
    imap(adjust_grid_axes, rows = 3, cols = 3, no_titles = "xy") %>%
    wrap_plots(ncol = 3, byrow = TRUE)

all_stoch_alt_layout <- "#24
                     135
                     #66"


all_stoch_alt_p_main <- wrap_elements(textGrob("Partitioning investment", rot = 90),
              clip = FALSE) +
    wrap_elements(textGrob("strong sub-additive", y = 1, vjust = 1), clip = FALSE) +
    all_sub_alt_p +
    wrap_elements(textGrob("strong super-additive", y = 1, vjust = 1), clip = FALSE) +
    all_super_alt_p +
    wrap_elements(textGrob("Conflict investment"), clip = FALSE) +
    plot_layout(heights = c(0.02, 1, 0.005),
                widths = c(0.01, 1, 1),
                design = all_stoch_alt_layout)


all_stoch_alt_p <- function() {
    grid.newpage()
    grid.draw(all_stoch_alt_p_main)
    grid.draw(textGrob("A", x = 0.01, y = 0.99, just = c("left", "top")))
    grid.draw(textGrob("B", x = 0.52, y = 0.99, just = c("center", "top")))
    invisible(NULL)
}

# all_stoch_alt_p()


# if (.RESAVE_PLOTS) save_plot(all_stoch_alt_p, 6.5, 3.7, .prefix = "S1-", .png = TRUE)