_results/fig5.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 <- FALSE
spp <- function(.n, .x = 0, .y = 1) {
rnds <- runif(.n, min = .x, max = .y)
M <- unname(rbind(rnds, rep(0, .n)))
M[,1:(.n %/% 2)] <- M[2:1,1:(.n %/% 2)]
M <- M[,sample.int(.n)]
return(M)
}
# Takes ~1 min
t0 <- Sys.time()
excl_sims <- mapply(.dp = c(0.1, 0.9), .big_seed = c(464714778, 1632530173),
FUN = function(.dp, .big_seed) {
set.seed(.big_seed)
ex_seeds <- sample.int(2^31-1, 1000)
ex_sims <- mclapply(1:1000, function(.r) {
set.seed(ex_seeds[.r])
tibble(eta = 0.6, d1 = -0.5, d2 = .dp,
final_t = 20e3, spp_gap_t = 5e3L,
save_every = 0L,
V0 = list(spp(20))) %>%
run_determ_1rep() %>%
select(c(spp, N, V1, V2)) %>%
unnest(c(spp, N, V1, V2)) %>%
mutate(rep = .r, seed = ex_seeds[.r])
}) %>%
do.call(what = bind_rows)
ex_sims$dp <- .dp
return(ex_sims)
},
SIMPLIFY = FALSE) %>%
do.call(what = bind_rows) %>%
select(dp, rep, seed, everything())
t1 <- Sys.time()
t1 - t0
# excl_sims
excl_p_df <- excl_sims %>%
group_by(dp, rep) %>%
summarize(v1 = sum(V1 > 0),
v2 = sum(V2 > 0),
n_spp = n(),
.groups = "drop") %>%
{group_by(., dp) %>% summarize(p1 = mean(v1 / n_spp)) %>% print(); .} %>%
{.[["v1"]] %>% table() %>% print(); .} %>%
{.[["n_spp"]] %>% table() %>% print(); .} %>%
# group_by(dp, n_spp) %>%
# summarize(across(c(v1), list(mean_p = ~ mean(.x / n_spp),
# mean_n = ~ mean(.x),
# min_n = ~ min(.x),
# max_n = ~ max(.x))),
# n = n())
group_by(dp, n_spp, v1) %>%
summarize(count = n(),
.groups = "drop") %>%
group_by(dp, n_spp) %>%
mutate(p = count / sum(count)) %>%
ungroup() %>%
mutate(p1 = v1 / n_spp, dp = factor(dp))
excl_p_df %>%
ggplot(aes(n_spp, p1)) +
geom_hline(yintercept = 0.5, linetype = "22", color = "gray70") +
# geom_bar(aes(fill = v1), position="stack", stat="identity") +
geom_point(aes(size = count)) +
# geom_point(data = excl_p_df %>%
# group_by(dp, n_spp) %>%
# summarize(p1 = sum(p1 * p), .groups = "drop"),
# color = "red") +
# stat_summary(fun = mean, geom = "line") +
xlab("Number of surviving species") +
ylab("Proportion investing in conflict traits") +
facet_wrap(~ dp)
exc_sims %>% filter(rep == 1)
set.seed(643327297)
tibble(eta = 0.6, d1 = -0.5, d2 = .dp,
final_t = 10e3, spp_gap_t = 5e3L,
save_every = 0L,
V0 = list(spp(10))) %>%
run_determ_1rep() %>%
select(c(spp, N, V1, V2)) %>%
unnest(c(spp, N, V1, V2))
# ----------------------------------------------------------------------------*
# ----------------------------------------------------------------------------*
# ----------------------------------------------------------------------------*
# ----------------------------------------------------------------------------*
# ----------------------------------------------------------------------------*
# ----------------------------------------------------------------------------*
#'
#' Simulate 2-species equilibrium community.
#'
one_dx_dp_equil_sim <- function(.dx, .dp, .final_t = 30e3L) {
# .dx = 0.6; .dp = 0.6
# rm(sim0, raw_sim0, .V0_0, .d_vec, .dx, .dp)
.d_vec <- c(- .dx, .dp)
.V0_0 <- rbind(cos(c(6, 84) * pi / 180) * 1,
sin(c(6, 84) * pi / 180) * 1) %>%
list()
raw_sim0 <- tibble(eta = 0, d1 = .d_vec[1], d2 = .d_vec[2],
final_t = .final_t, spp_gap_t = 0L, save_every = 1L,
V0 = .V0_0) %>%
run_determ_1rep()
# Stable equilibrium?
if (!equil_check(raw_sim0)) {
stop("\nDidn't reach equilibrium for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
}
if ((max(Re(raw_sim0[["eig"]][[1]])) - 1) > 1e-6) {
stop("\nEquilibrium not (neutrally) stable for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
}
sim0 <- raw_sim0 %>%
select(spp, time, N, V1, V2) %>%
unnest(c(spp, time, N, V1, V2)) %>%
filter(time == max(time)) %>%
mutate(Omega = calc_omega(., d_vec = .d_vec)) %>%
select(-time) %>%
mutate(invest = sqrt(V1^2 + V2^2))
if (nrow(sim0) != ncol(.V0_0[[1]])) {
stop("\nNot all species survived initial sim for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
}
out <- tibble(dx = .dx, dp = .dp, sim0 = list(sim0))
return(out)
}
#' Simulations for effects on conflicter's investment.
#' Takes ~90 sec
# invest_d_sims <- crossing(.dx = seq(0.1, 1.5, 0.1) %>% round(1),
# .dp = seq(0.1, 1.5, 0.05) %>% round(2)) %>%
# as.list() %>%
# c(list(FUN = one_dx_dp_equil_sim, SIMPLIFY = FALSE)) %>%
# do.call(what = mcmapply) %>%
# do.call(what = bind_rows) %>%
# mutate(invest_c = map_dbl(sim0, ~ .x$invest[[1]]),
# invest_p = map_dbl(sim0, ~ .x$invest[[2]]))
# saveRDS(invest_d_sims, "_results/results-rds/invest_d_sims.rds")
invest_d_sims <- readRDS("_results/results-rds/invest_d_sims.rds")
#'
#' How does investment by the conflict / partitioning investor change
#' as we change dp and dx?
#'
dp_invest_p <- invest_d_sims %>%
filter(dx %in% c(0.3, 0.5, 0.7, 0.9)) %>%
filter(dp <= 0.9) %>%
ggplot(aes(dp, invest_c, color = factor(dx))) +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_line(size = 0.5) +
geom_point(size = 1) +
geom_line(aes(y = invest_p), size = 0.5, linetype = "21") +
geom_point(aes(y = invest_p), size = 1, shape = 15) +
geom_text(data = tibble(dp = c(0.1, 0.85),
invest_c = c(0.55, 1.1),
lab = c("conflicter", "partitioner")),
aes(label = lab), hjust = c(0, 1), vjust = 1, color = "black") +
scale_color_manual(expression(d[x]),
values = c("gray80", "gray60", "gray40", "black")) +
scale_x_continuous(expression(d[p]), breaks = c(1,5,9)/10) +
scale_y_continuous("Equilibrium investment", breaks = 0.5 * 0:2,
limits = c(0, 1.35)) +
NULL
dp_omega_p <- invest_d_sims %>%
filter(dx %in% c(0.3, 0.5, 0.7, 0.9)) %>%
filter(dp <= 0.9) %>%
mutate(omega_c = map_dbl(sim0, ~ .x$Omega[[1]]),
omega_p = map_dbl(sim0, ~ .x$Omega[[2]])) %>%
ggplot(aes(dp, omega_c, color = factor(dx))) +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_line(size = 0.5) +
geom_point(size = 1) +
geom_line(aes(y = omega_p), size = 0.5, linetype = "21") +
geom_point(aes(y = omega_p), size = 1, shape = 15) +
scale_color_manual(expression(d[x]),
values = c("gray80", "gray60", "gray40", "black")) +
scale_x_continuous(expression(d[p]), breaks = c(1,5,9)/10) +
scale_y_continuous("Equilibrium ECN", breaks = 2500 * 0:2,
limits = c(0, 6000)) +
NULL
dp_abund_p <- invest_d_sims %>%
filter(dx %in% c(0.3, 0.5, 0.7, 0.9)) %>%
filter(dp <= 0.9) %>%
mutate(N_c = map_dbl(sim0, ~ .x$N[[1]]),
N_p = map_dbl(sim0, ~ .x$N[[2]])) %>%
ggplot(aes(dp, N_c, color = factor(dx))) +
geom_line(size = 0.5) +
geom_point(size = 1) +
geom_line(aes(y = N_p), size = 0.5, linetype = "21") +
geom_point(aes(y = N_p), size = 1, shape = 15) +
scale_color_manual(expression(d[x]),
values = c("gray80", "gray60", "gray40", "black")) +
scale_x_continuous(expression(d[p]), breaks = c(1,5,9)/10) +
scale_y_continuous("Equilibrium abundance", breaks = 1000 * 2:4,
limits = c(1900, 4305)) +
NULL
dp_invest_abund_p <- dp_invest_p + dp_omega_p + dp_abund_p +
plot_annotation(tag_levels = "A") +
plot_layout(ncol = 1, guides = "collect")
# if (.RESAVE_PLOTS) save_plot(dp_invest_abund_p, 4, 8, .prefix = "S4-")
#' These are similar to plots above, except that dx == 0.5 for simplicity,
#' plus they're simpler overall bc they'll be put together with others
#' to make Figure 5 using Illustrator.
#'
invest_d_p <- invest_d_sims %>%
filter(dx == 0.5) %>%
filter(dp <= 0.9) %>%
select(-sim0) %>%
pivot_longer(starts_with("invest"), names_to = "spp",
values_to = "invest") %>%
mutate(spp = ifelse(str_ends(spp, "_c"), "conflicter", "partitioner") %>%
factor(levels = c("conflicter", "partitioner"))) %>%
ggplot(aes(dp, invest, color = spp)) +
# geom_vline(xintercept = c(0.3, 0.7), size = 0.5, color = "gray70",
# linetype = "22") +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_line(size = 0.5) +
geom_point(size = 1.5) +
scale_color_manual(values = plasma(100)[c(1, 85, 65)][-3]) +
scale_x_continuous(expression(d[p]), breaks = c(1,3,5,7,9)/10) +
scale_y_continuous("Equilibrium investment", breaks = 0.5 * 0:2,
limits = c(0, 1.35)) +
theme(axis.title = element_blank(),
axis.text = element_blank(),
legend.position = "none",
plot.margin = margin(0,0,0,0)) +
NULL
if (.RESAVE_PLOTS) save_plot(invest_d_p, 1.75, 1,
.name = "fig5_sub-plots/invest_d_p")
abund_d_p <- invest_d_sims %>%
filter(dx == 0.5) %>%
filter(dp <= 0.9) %>%
mutate(N_c = map_dbl(sim0, ~ .x$N[[1]]),
N_p = map_dbl(sim0, ~ .x$N[[2]])) %>%
select(-sim0, -starts_with("invest")) %>%
pivot_longer(starts_with("N"), names_to = "spp",
values_to = "N") %>%
mutate(spp = ifelse(str_ends(spp, "_c"), "conflicter", "partitioner") %>%
factor(levels = c("conflicter", "partitioner"))) %>%
ggplot(aes(dp, N, color = spp)) +
# geom_vline(xintercept = c(0.3, 0.7), size = 0.5, color = "gray70",
# linetype = "22") +
geom_line(size = 0.5) +
geom_point(size = 1.5) +
scale_color_manual(values = plasma(100)[c(1, 85, 65)][-3]) +
scale_x_continuous(expression(d[p]), breaks = c(1,3,5,7,9)/10) +
scale_y_continuous("Equilibrium abundance", breaks = 2400 + 800 * 0:2,
limits = c(2400, 4305)) +
theme(axis.title = element_blank(),
axis.text = element_blank(),
legend.position = "none",
plot.margin = margin(0,0,0,0)) +
NULL
if (.RESAVE_PLOTS) save_plot(abund_d_p, 1.75, 1,
.name = "fig5_sub-plots/abund_d_p")
#'
#' Simulate invaders into a 2-species community.
#'
one_dx_dp_inv_sim <- function(.dx, .dp, .sim0, .sim_max_t) {
# .dx = 0.6; .dp = 0.6; .sim_max_t = 10e3L
# rm(out, survs, x_vals, O_inv, sim0, raw_sim0, .V0_0, .d_vec,
# .dx, .dp, .sim_max_t)
stopifnot(!inherits(.sim0, "list"))
.d_vec <- c(- .dx, .dp)
O_inv <- with(.sim0, sum(N * exp(- .d_vec[1] * V1^2 - .d_vec[2] * V2^2)))
x_vals <- seq(0.1, 2.5, 0.1) %>%
round(1)
survs <- x_vals %>%
map(function(.x) {
#' No point in doing simulation if we know the invader's fitness
#' will be < 1
f0 <- with(formals(quant_gen), {
exp(r0 - f * .x^2 - a0 - a0 * exp(-.x^2) * O_inv)
})
if (f0 < 1) return(c(rep(1L, nrow(.sim0)), 0L))
.N0 <- c(.sim0$N, 1)
.V0 <- cbind(rbind(.sim0$V1, .sim0$V2), c(.x, 0)) %>%
list()
raw_sim1 <- tibble(eta = 0, d1 = .d_vec[1], d2 = .d_vec[2],
final_t = .sim_max_t, spp_gap_t = 0L,
save_every = 1L,
V0 = .V0, N0 = list(.N0)) %>%
run_determ_1rep()
surv_spp <- raw_sim1 %>%
select(spp, time, N, V1, V2) %>%
unnest(c(spp, time, N, V1, V2)) %>%
filter(time == max(time)) %>%
.[["spp"]] %>%
as.integer()
surv_vec <- as.integer(1:(nrow(.sim0)+1) %in% surv_spp)
# Stable equilibrium?
if (!equil_check(raw_sim1)) {
warning("\nDidn't reach equilibrium at `x` = ", .x,
" for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
surv_vec[surv_vec == 1] <- NA_integer_
}
if ((max(Re(raw_sim1[["eig"]][[1]])) - 1) > 1e-6) {
warning("\nEquilibrium unstable at `x` = ", .x,
" for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
surv_vec[surv_vec == 1] <- NA_integer_
}
return(surv_vec)
})
out <- cbind(x_vals, do.call(rbind, survs)) %>%
as.data.frame() %>%
set_names(c("x", "spp1", "spp2", "spp3")) %>%
as_tibble() %>%
mutate(dx = .dx, dp = .dp) %>%
select(dx, dp, everything())
return(out)
}
# # Takes ~9 min
# dx_dp_sims <- invest_d_sims %>%
# select(dx, dp, sim0) %>%
# rename_with(~ paste0(".", .x)) %>%
# # Some combos need to be run longer to reach equilibrium:
# mutate(.sim_max_t = case_when(.dp >= 0.7 & .dp <= 0.85 & .dx >= 0.7 ~ 30e3L,
# .dp >= 0.7 & .dp <= 0.85 & .dx > 0.6 ~ 20e3L,
# .dp >= 0.7 & .dp <= 0.85 & .dx >= 0.2 ~ 30e3L,
# TRUE ~ 10e3L)) %>%
# # head() %>%
# as.list() %>%
# c(list(FUN = one_dx_dp_inv_sim, SIMPLIFY = FALSE)) %>%
# do.call(what = mcmapply) %>%
# do.call(what = bind_rows)
# saveRDS(dx_dp_sims, "_results/results-rds/dx_dp_sims.rds")
dx_dp_sims <- readRDS("_results/results-rds/dx_dp_sims.rds")
#' This should return no rows bc it checks for any reps not reaching an
#' equilibrium or reach an unstable equilibrium
dx_dp_sims %>%
filter(if_any(starts_with("spp"), is.na))
#' Only 4 possible outcomes:
#' 1. both residents persist, invader excluded
#' 2. species 1 excluded, invader and species 2 coexist
#' 3. all 3 species persist
#' 4. both residents excluded, invader persists
dx_dp_sims %>%
mutate(outcome = paste0(spp1, spp2, spp3)) %>%
distinct(outcome)
dx_dp_sims %>% filter(dx == 0.5 & (dp == 0.3 | dp == 0.7) & x == 1)
invest_d_comm_p <- invest_d_sims %>%
filter(dx == 0.5, dp %in% c(0.3, 0.7)) %>%
arrange(dp) %>%
.[["sim0"]] %>%
map(function(dd) {
dd %>%
ggplot(aes(V1, V2, color = spp))+
geom_abline(slope = 1, intercept = 0,
linetype = 2, color = "gray70") +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_vline(xintercept = 0, color = "gray70", size = 0.5) +
geom_point(shape = 19, size = 3) +
scale_color_manual(values = plasma(100)[c(1, 85)], guide = "none") +
scale_y_continuous("Partitioning investment", breaks = 0.5 * 0:2) +
scale_x_continuous("Conflict investment", breaks = 0.5 * 0:2) +
coord_equal(xlim = c(-0.05, 1.3), ylim = c(-0.05, 1.3),
clip = "off") +
theme(panel.background = element_blank(),
axis.title = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
plot.margin = margin(0,0,0,0)) +
NULL
})
if (.RESAVE_PLOTS) {
save_plot(invest_d_comm_p[[1]], 0.8, 0.8, .name = "fig5_sub-plots/comm_0.3")
save_plot(invest_d_comm_p[[2]], 0.8, 0.8, .name = "fig5_sub-plots/comm_0.7")
}
invade_dp <- function(.dp, .inv_x = 1) {
# .dp = 0.7; .inv_x = 1
# rm(.dx, .dp, .inv_x, equil, raw_sim1, prequel, out)
.dx <- 0.5
stopifnot(length(.dp) == 1 && is.numeric(.dp) && .dp > 0)
stopifnot(length(.inv_x) == 1 && is.numeric(.inv_x) && .inv_x > 0)
equil <- invest_d_sims %>%
filter(dx == .dx, dp == .dp) %>%
.[["sim0"]] %>%
.[[1]]
raw_sim1 <- tibble(eta = 0, d1 = -.dx, d2 = .dp,
final_t = 10e3L, spp_gap_t = 0L,
save_every = 1L,
V0 = list(rbind(c(equil$V1, .inv_x),
c(equil$V2, 0))),
N0 = list(c(equil$N, 1))) %>%
run_determ_1rep() %>%
select(spp, time, N, V1, V2) %>%
unnest(c(spp, time, N, V1, V2)) %>%
mutate(Omega = calc_omega(., d_vec = c(-.dx, .dp)),
invest = sqrt(V1^2 + V2^2),
dp = .dp)
prequel <- equil %>% mutate(dp = .dp, time = -1L)
out <- bind_rows(prequel, raw_sim1) %>%
select(dp, everything())
return(out)
}
invade_dp_df <- map_dfr(c(0.3, 0.7), invade_dp)
invade_dp_p_list <- sort(unique(invade_dp_df$dp)) %>%
map(function(.dp) {
.dd <- invade_dp_df %>% filter(dp == .dp)
ext_df <- .dd %>%
group_by(spp) %>%
filter(time == max(time)) %>%
ungroup() %>%
filter(time < max(time))
y_scales <- list(invest = scale_y_continuous("Total investment",
limits = c(-0.16, 1.6),
breaks = c(0, 0.7, 1.4)),
Omega = scale_y_continuous("ECN",
limits = c(-1100, 11000),
breaks = 5e3*0:2),
N = scale_y_continuous("Abundance",
limits = c(-410, 4100),
breaks = 2e3*0:2))
plots <- imap(y_scales, function(y_scale, y_var) {
.dd %>%
filter(time < 350) %>%
mutate(time = ifelse(time < 0, time * 100, time)) %>%
ggplot(aes(time, !!sym(y_var), color = spp)) +
geom_vline(xintercept = 0, color = "gray70", size = 0.5) +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_line(size = 0.75) +
geom_point(data = ext_df, size = 2, shape = 4) +
scale_color_manual(values = plasma(100)[c(1, 85, 65)],
guide = "none") +
scale_x_continuous("Time", breaks = 100 * 0:3) +
y_scale +
theme(axis.title = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_blank(),
plot.margin = margin(0,0,0,0)) +
coord_cartesian(xlim = c(-50, NA), expand = FALSE)
})
p <- wrap_plots(plots, ncol = 1) &
theme(plot.margin = margin(0,0,0,0))
return(p)
})
# wrap_plots(invade_dp_p_list, ncol = 2)
if (.RESAVE_PLOTS) {
save_plot(invade_dp_p_list[[1]], 2, 1.75, .name = "fig5_sub-plots/invade_0.3")
save_plot(invade_dp_p_list[[2]], 2, 1.75, .name = "fig5_sub-plots/invade_0.7")
}
# =============================================================================*
# =============================================================================*
# Full outcomes
# =============================================================================*
# =============================================================================*
d_outcome_p <- dx_dp_sims %>%
mutate(outcome = case_when(spp1 == 1 & spp3 == 1 ~ "3-spp coexist",
spp1 == 1 & spp3 == 0 ~ "invader excluded",
spp1 == 0 & spp2 == 1 & spp3 == 1 ~ "conflict resident excluded",
spp1 == 0 & spp2 == 0 & spp3 == 1 ~ "both residents excluded",
TRUE ~ NA_character_) %>%
factor(levels = c("both residents excluded",
"conflict resident excluded",
"3-spp coexist",
"invader excluded"))) %>%
mutate(dx_f = factor(dx, levels = sort(unique(dx)),
labels = paste("d[x] ==", sort(unique(dx))))) %>%
ggplot(aes(dp, x)) +
geom_raster(aes(fill = outcome)) +
scale_fill_manual(NULL, values = c("gray80", plasma(100)[c(70, 40, 20)])) +
scale_y_continuous(expression("Invader" ~ x[0]), breaks = 0:2) +
scale_x_continuous(expression(d[p]), breaks = c(2,8,14)/10) +
facet_wrap(~ dx_f, labeller = label_parsed, nrow = 3) +
guides(fill = guide_legend(ncol = 1)) +
theme(legend.position = "right") +
NULL
if (.RESAVE_PLOTS) save_plot(d_outcome_p, 6, 4, .prefix = "S5-", .png = TRUE)
dx_dp_minx <- dx_dp_sims %>%
split(interaction(.$dx, .$dp)) %>%
map_dfr(function(.d) {
.res_excl <- .d %>%
filter(spp1 == 0 & spp3 == 1)
if (nrow(.res_excl) == 0) {
.min_x <- 100
} else {
.min_x <- min(.res_excl$x)
}
out <- select(.d[1,], dx, dp) %>%
mutate(min_x = .min_x)
return(out)
}) %>%
mutate(inf = factor(min_x == 100),
min_x = ifelse(min_x == 100, 2.2, min_x))
dp_minx_p <- dx_dp_minx %>%
filter(dx %in% (c(3, 5, 7, 9) / 10)) %>%
filter(dp <= 0.9) %>%
group_by(dx) %>%
filter(inf == "FALSE" | dp ==
if (any(inf == "TRUE")) max(dp[inf == "TRUE"]) else -Inf) %>%
ungroup() %>%
ggplot(aes(dp, min_x)) +
geom_hline(yintercept = 0.1, color = "gray70", size = 0.5) +
geom_hline(yintercept = 2.2, color = "gray70", size = 0.5) +
geom_line(aes(color = factor(dx)), size = 0.5) +
geom_point(aes(shape = inf, color = factor(dx)), size = 1.5) +
scale_shape_manual(values = c(19, 1), guide = "none") +
# scale_color_viridis_d(expression(d[x]), begin = 0.2, end = 0.8) +
scale_color_manual(values = c("gray80", "gray60", "gray40", "black")) +
scale_x_continuous(expression(d[p]), breaks = c(1,5,9)/10) +
scale_y_continuous(expression(x[inv]), breaks = 0:2) +
theme(axis.title = element_blank(),
axis.text = element_blank(),
legend.position = "none",
plot.margin = margin(0,0,0,0)) +
NULL
if (.RESAVE_PLOTS) save_plot(dp_minx_p, 2, 1.75,
.name = "fig5_sub-plots/dp_minx")
paired_invasions <- function(.dx, .dp, .steps = c(-0.15, 0.15)) {
# .dx = 0.7; .dp = 0.6; .step = 0.15
# rm(.d_vec, .min_x, .x_vec, equil, out_list, out, i, .x, raw_sim1, prequel)
stopifnot(length(.steps) == 2 && is.numeric(.steps) &&
.steps[1] < 0 && .steps[2] > 0)
.d_vec = c(- abs(.dx), abs(.dp))
.min_x <- dx_dp_minx %>%
filter(dx == abs(.d_vec[1]), dp == .d_vec[2]) %>%
.[["min_x"]]
.x_vec <- .min_x + .steps
equil <- invest_d_sims %>%
filter(dx == abs(.d_vec[1]), dp == .d_vec[2]) %>%
.[["sim0"]] %>%
.[[1]]
out_list <- rep(list(NA), 2)
for (i in 1:length(.x_vec)) {
.x = .x_vec[i]
raw_sim1 <- tibble(eta = 0, d1 = .d_vec[1], d2 = .d_vec[2],
final_t = 10e3L, spp_gap_t = 0L,
save_every = 1L,
V0 = list(rbind(c(equil$V1, .x),
c(equil$V2, 0))),
N0 = list(c(equil$N, 1))) %>%
run_determ_1rep()
# Stable equilibrium?
if (!equil_check(raw_sim1)) {
stop("\nDidn't reach equilibrium at `x` = ", .x,
" for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
}
if ((max(Re(raw_sim1[["eig"]][[1]])) - 1) > 1e-6) {
stop("\nEquilibrium unstable at `x` = ", .x,
" for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
}
out <- raw_sim1 %>%
select(spp, time, N, V1, V2) %>%
unnest(c(spp, time, N, V1, V2)) %>%
mutate(Omega = calc_omega(., d_vec = .d_vec),
invest = sqrt(V1^2 + V2^2),
x_inv = .x) %>%
select(x_inv, everything())
prequel <- equil %>% mutate(x_inv = .x, time = -1L)
out_list[[i]] <- bind_rows(prequel, out)
}
do.call(bind_rows, out_list)
}
# The below code shows plots of the weird edge case issue in Figure S5.
# x <- paired_invasions(0.1, 1.35, .steps = c(-2, 0.5)) %>%
# filter(x_inv == min(x_inv))
# y <- paired_invasions(0.3, 1.35, .steps = c(-0.01, 0.4)) %>%
# filter(x_inv == max(x_inv))
# # The problem one:
# z <- paired_invasions(0.2, 1.35, .steps = c(-0.6, 0.9)) %>%
# filter(x_inv == min(x_inv))
#
# problem_plot <- function(.d, .dx, y_var, .max, .leg = TRUE) {
# p <- .d %>%
# ggplot(aes(time, !!sym(y_var), color = spp)) +
# ggtitle(sprintf("dx = %.2f", .dx)) +
# geom_line() +
# scale_y_continuous(limits = c(0, .max)) +
# scale_color_manual(values = plasma(100)[c(1, 85, 65)])
# if (!.leg) p <- p +
# theme(legend.position = "none")
# return(p)
# }
#
#
# problem_plot(x, 0.1, "invest", 1.4) +
# problem_plot(z, 0.2, "invest", 1.4) +
# problem_plot(y, 0.3, "invest", 1.4) +
# problem_plot(x, 0.1, "Omega", 7000) +
# problem_plot(z, 0.2, "Omega", 7000) +
# problem_plot(y, 0.3, "Omega", 7000) +
# problem_plot(x, 0.1, "N", 4700) +
# problem_plot(z, 0.2, "N", 4700) +
# problem_plot(y, 0.3, "N", 4700) +
# plot_layout(ncol = 3, byrow = TRUE, guides = "collect")
paired_inv_df <- paired_invasions(0.5, 0.6, .steps = c(-0.9, 0.1))
paired_inv_p <- sort(unique(paired_inv_df$x_inv)) %>%
map(function(xinv) {
.dd <- paired_inv_df %>% filter(x_inv == xinv)
ext_df <- .dd %>%
group_by(spp) %>%
filter(time == max(time)) %>%
ungroup() %>%
filter(time < max(time))
y_scales <- list(invest = scale_y_continuous("Total investment",
limits = c(-0.17, 1.7),
breaks = c(0, 0.8, 1.6)),
Omega = scale_y_continuous("ECN",
limits = c(-900, 9000),
breaks = 4e3*0:2),
N = scale_y_continuous("Abundance",
limits = c(-410, 4100),
breaks = 2e3*0:2))
plots <- imap(y_scales, function(y_scale, y_var) {
.dd %>%
filter(time < 400) %>%
mutate(time = ifelse(time < 0, time * 100, time)) %>%
ggplot(aes(time, !!sym(y_var), color = spp)) +
geom_vline(xintercept = 0, color = "gray70", size = 0.5) +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_line(size = 0.75) +
geom_point(data = ext_df, size = 2, shape = 4) +
scale_color_manual(values = plasma(100)[c(1, 85, 65)],
guide = "none") +
scale_x_continuous("Time", breaks = 100 * 0:3) +
y_scale +
theme(axis.title = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_blank(),
plot.margin = margin(0,0,0,0)) +
coord_cartesian(xlim = c(-50, NA), expand = FALSE)
})
p <- wrap_plots(plots, ncol = 1) &
theme(plot.margin = margin(0,0,0,0))
return(p)
})
# wrap_plots(paired_inv_p, ncol = 2)
if (.RESAVE_PLOTS) {
save_plot(paired_inv_p[[1]], 2, 1.75, .name = "fig5_sub-plots/under")
save_plot(paired_inv_p[[2]], 2, 1.75, .name = "fig5_sub-plots/over")
}
paired_inv_comm_p <- invest_d_sims %>%
filter(dx == 0.5, dp == 0.6) %>%
.[["sim0"]] %>%
.[[1]] %>%
ggplot(aes(V1, V2, color = spp))+
geom_abline(slope = 1, intercept = 0,
linetype = 2, color = "gray70") +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_vline(xintercept = 0, color = "gray70", size = 0.5) +
geom_point(shape = 19, size = 3) +
# geom_point(shape = 19, size = 3, alpha = 0.75) +
# geom_point(shape = 1, size = 3) +
geom_segment(y = -0.05, yend = 0.05,
x = dx_dp_minx %>%
filter(dx == 0.5, dp == 0.6) %>%
.[["min_x"]],
xend = dx_dp_minx %>%
filter(dx == 0.5, dp == 0.6) %>%
.[["min_x"]], color = "black", size = 0.5) +
scale_color_manual(values = plasma(100)[c(1, 85)], guide = "none") +
scale_y_continuous("Partitioning investment", breaks = 0:2) +
scale_x_continuous("Conflict investment", breaks = 0:2) +
coord_equal(xlim = c(-0.025, 2), ylim = c(-0.025, 2), clip = "off") +
theme(panel.background = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
plot.margin = margin(0,0,0,0))
paired_inv_comm_p_invasion <-
list("min", "max") %>%
set_names() %>%
map(~ eval(parse(text = .x))) %>%
map(function(.f) {
.dd <- filter(paired_inv_df, time == 0, spp == 3) %>%
filter(V1 == .f(V1))
p <- paired_inv_comm_p +
geom_point(data = .dd,
shape = 19, size = 3, color = plasma(100)[c(65)])
gt <- ggplot_gtable(ggplot_build(p))
gt$layout$z[grepl("^panel", gt$layout$name)] <- max(gt$layout$z) + 1
gt
})
if (.RESAVE_PLOTS) {
save_plot(paired_inv_comm_p, 1, 1, .name = "fig5_sub-plots/comm")
save_plot(function() {grid.newpage(); grid.draw(paired_inv_comm_p_invasion$max)},
0.8, 0.8, .name = "fig5_sub-plots/comm-above")
save_plot(function() {grid.newpage(); grid.draw(paired_inv_comm_p_invasion$min)},
0.8, 0.8, .name = "fig5_sub-plots/comm-below")
}
lucasnell/evo_alt_states documentation built on Aug. 17, 2022, 5:34 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# 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 <- FALSE
spp <- function(.n, .x = 0, .y = 1) {
rnds <- runif(.n, min = .x, max = .y)
M <- unname(rbind(rnds, rep(0, .n)))
M[,1:(.n %/% 2)] <- M[2:1,1:(.n %/% 2)]
M <- M[,sample.int(.n)]
return(M)
}
# Takes ~1 min
t0 <- Sys.time()
excl_sims <- mapply(.dp = c(0.1, 0.9), .big_seed = c(464714778, 1632530173),
FUN = function(.dp, .big_seed) {
set.seed(.big_seed)
ex_seeds <- sample.int(2^31-1, 1000)
ex_sims <- mclapply(1:1000, function(.r) {
set.seed(ex_seeds[.r])
tibble(eta = 0.6, d1 = -0.5, d2 = .dp,
final_t = 20e3, spp_gap_t = 5e3L,
save_every = 0L,
V0 = list(spp(20))) %>%
run_determ_1rep() %>%
select(c(spp, N, V1, V2)) %>%
unnest(c(spp, N, V1, V2)) %>%
mutate(rep = .r, seed = ex_seeds[.r])
}) %>%
do.call(what = bind_rows)
ex_sims$dp <- .dp
return(ex_sims)
},
SIMPLIFY = FALSE) %>%
do.call(what = bind_rows) %>%
select(dp, rep, seed, everything())
t1 <- Sys.time()
t1 - t0
# excl_sims
excl_p_df <- excl_sims %>%
group_by(dp, rep) %>%
summarize(v1 = sum(V1 > 0),
v2 = sum(V2 > 0),
n_spp = n(),
.groups = "drop") %>%
{group_by(., dp) %>% summarize(p1 = mean(v1 / n_spp)) %>% print(); .} %>%
{.[["v1"]] %>% table() %>% print(); .} %>%
{.[["n_spp"]] %>% table() %>% print(); .} %>%
# group_by(dp, n_spp) %>%
# summarize(across(c(v1), list(mean_p = ~ mean(.x / n_spp),
# mean_n = ~ mean(.x),
# min_n = ~ min(.x),
# max_n = ~ max(.x))),
# n = n())
group_by(dp, n_spp, v1) %>%
summarize(count = n(),
.groups = "drop") %>%
group_by(dp, n_spp) %>%
mutate(p = count / sum(count)) %>%
ungroup() %>%
mutate(p1 = v1 / n_spp, dp = factor(dp))
excl_p_df %>%
ggplot(aes(n_spp, p1)) +
geom_hline(yintercept = 0.5, linetype = "22", color = "gray70") +
# geom_bar(aes(fill = v1), position="stack", stat="identity") +
geom_point(aes(size = count)) +
# geom_point(data = excl_p_df %>%
# group_by(dp, n_spp) %>%
# summarize(p1 = sum(p1 * p), .groups = "drop"),
# color = "red") +
# stat_summary(fun = mean, geom = "line") +
xlab("Number of surviving species") +
ylab("Proportion investing in conflict traits") +
facet_wrap(~ dp)
exc_sims %>% filter(rep == 1)
set.seed(643327297)
tibble(eta = 0.6, d1 = -0.5, d2 = .dp,
final_t = 10e3, spp_gap_t = 5e3L,
save_every = 0L,
V0 = list(spp(10))) %>%
run_determ_1rep() %>%
select(c(spp, N, V1, V2)) %>%
unnest(c(spp, N, V1, V2))
# ----------------------------------------------------------------------------*
# ----------------------------------------------------------------------------*
# ----------------------------------------------------------------------------*
# ----------------------------------------------------------------------------*
# ----------------------------------------------------------------------------*
# ----------------------------------------------------------------------------*
#'
#' Simulate 2-species equilibrium community.
#'
one_dx_dp_equil_sim <- function(.dx, .dp, .final_t = 30e3L) {
# .dx = 0.6; .dp = 0.6
# rm(sim0, raw_sim0, .V0_0, .d_vec, .dx, .dp)
.d_vec <- c(- .dx, .dp)
.V0_0 <- rbind(cos(c(6, 84) * pi / 180) * 1,
sin(c(6, 84) * pi / 180) * 1) %>%
list()
raw_sim0 <- tibble(eta = 0, d1 = .d_vec[1], d2 = .d_vec[2],
final_t = .final_t, spp_gap_t = 0L, save_every = 1L,
V0 = .V0_0) %>%
run_determ_1rep()
# Stable equilibrium?
if (!equil_check(raw_sim0)) {
stop("\nDidn't reach equilibrium for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
}
if ((max(Re(raw_sim0[["eig"]][[1]])) - 1) > 1e-6) {
stop("\nEquilibrium not (neutrally) stable for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
}
sim0 <- raw_sim0 %>%
select(spp, time, N, V1, V2) %>%
unnest(c(spp, time, N, V1, V2)) %>%
filter(time == max(time)) %>%
mutate(Omega = calc_omega(., d_vec = .d_vec)) %>%
select(-time) %>%
mutate(invest = sqrt(V1^2 + V2^2))
if (nrow(sim0) != ncol(.V0_0[[1]])) {
stop("\nNot all species survived initial sim for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
}
out <- tibble(dx = .dx, dp = .dp, sim0 = list(sim0))
return(out)
}
#' Simulations for effects on conflicter's investment.
#' Takes ~90 sec
# invest_d_sims <- crossing(.dx = seq(0.1, 1.5, 0.1) %>% round(1),
# .dp = seq(0.1, 1.5, 0.05) %>% round(2)) %>%
# as.list() %>%
# c(list(FUN = one_dx_dp_equil_sim, SIMPLIFY = FALSE)) %>%
# do.call(what = mcmapply) %>%
# do.call(what = bind_rows) %>%
# mutate(invest_c = map_dbl(sim0, ~ .x$invest[[1]]),
# invest_p = map_dbl(sim0, ~ .x$invest[[2]]))
# saveRDS(invest_d_sims, "_results/results-rds/invest_d_sims.rds")
invest_d_sims <- readRDS("_results/results-rds/invest_d_sims.rds")
#'
#' How does investment by the conflict / partitioning investor change
#' as we change dp and dx?
#'
dp_invest_p <- invest_d_sims %>%
filter(dx %in% c(0.3, 0.5, 0.7, 0.9)) %>%
filter(dp <= 0.9) %>%
ggplot(aes(dp, invest_c, color = factor(dx))) +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_line(size = 0.5) +
geom_point(size = 1) +
geom_line(aes(y = invest_p), size = 0.5, linetype = "21") +
geom_point(aes(y = invest_p), size = 1, shape = 15) +
geom_text(data = tibble(dp = c(0.1, 0.85),
invest_c = c(0.55, 1.1),
lab = c("conflicter", "partitioner")),
aes(label = lab), hjust = c(0, 1), vjust = 1, color = "black") +
scale_color_manual(expression(d[x]),
values = c("gray80", "gray60", "gray40", "black")) +
scale_x_continuous(expression(d[p]), breaks = c(1,5,9)/10) +
scale_y_continuous("Equilibrium investment", breaks = 0.5 * 0:2,
limits = c(0, 1.35)) +
NULL
dp_omega_p <- invest_d_sims %>%
filter(dx %in% c(0.3, 0.5, 0.7, 0.9)) %>%
filter(dp <= 0.9) %>%
mutate(omega_c = map_dbl(sim0, ~ .x$Omega[[1]]),
omega_p = map_dbl(sim0, ~ .x$Omega[[2]])) %>%
ggplot(aes(dp, omega_c, color = factor(dx))) +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_line(size = 0.5) +
geom_point(size = 1) +
geom_line(aes(y = omega_p), size = 0.5, linetype = "21") +
geom_point(aes(y = omega_p), size = 1, shape = 15) +
scale_color_manual(expression(d[x]),
values = c("gray80", "gray60", "gray40", "black")) +
scale_x_continuous(expression(d[p]), breaks = c(1,5,9)/10) +
scale_y_continuous("Equilibrium ECN", breaks = 2500 * 0:2,
limits = c(0, 6000)) +
NULL
dp_abund_p <- invest_d_sims %>%
filter(dx %in% c(0.3, 0.5, 0.7, 0.9)) %>%
filter(dp <= 0.9) %>%
mutate(N_c = map_dbl(sim0, ~ .x$N[[1]]),
N_p = map_dbl(sim0, ~ .x$N[[2]])) %>%
ggplot(aes(dp, N_c, color = factor(dx))) +
geom_line(size = 0.5) +
geom_point(size = 1) +
geom_line(aes(y = N_p), size = 0.5, linetype = "21") +
geom_point(aes(y = N_p), size = 1, shape = 15) +
scale_color_manual(expression(d[x]),
values = c("gray80", "gray60", "gray40", "black")) +
scale_x_continuous(expression(d[p]), breaks = c(1,5,9)/10) +
scale_y_continuous("Equilibrium abundance", breaks = 1000 * 2:4,
limits = c(1900, 4305)) +
NULL
dp_invest_abund_p <- dp_invest_p + dp_omega_p + dp_abund_p +
plot_annotation(tag_levels = "A") +
plot_layout(ncol = 1, guides = "collect")
# if (.RESAVE_PLOTS) save_plot(dp_invest_abund_p, 4, 8, .prefix = "S4-")
#' These are similar to plots above, except that dx == 0.5 for simplicity,
#' plus they're simpler overall bc they'll be put together with others
#' to make Figure 5 using Illustrator.
#'
invest_d_p <- invest_d_sims %>%
filter(dx == 0.5) %>%
filter(dp <= 0.9) %>%
select(-sim0) %>%
pivot_longer(starts_with("invest"), names_to = "spp",
values_to = "invest") %>%
mutate(spp = ifelse(str_ends(spp, "_c"), "conflicter", "partitioner") %>%
factor(levels = c("conflicter", "partitioner"))) %>%
ggplot(aes(dp, invest, color = spp)) +
# geom_vline(xintercept = c(0.3, 0.7), size = 0.5, color = "gray70",
# linetype = "22") +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_line(size = 0.5) +
geom_point(size = 1.5) +
scale_color_manual(values = plasma(100)[c(1, 85, 65)][-3]) +
scale_x_continuous(expression(d[p]), breaks = c(1,3,5,7,9)/10) +
scale_y_continuous("Equilibrium investment", breaks = 0.5 * 0:2,
limits = c(0, 1.35)) +
theme(axis.title = element_blank(),
axis.text = element_blank(),
legend.position = "none",
plot.margin = margin(0,0,0,0)) +
NULL
if (.RESAVE_PLOTS) save_plot(invest_d_p, 1.75, 1,
.name = "fig5_sub-plots/invest_d_p")
abund_d_p <- invest_d_sims %>%
filter(dx == 0.5) %>%
filter(dp <= 0.9) %>%
mutate(N_c = map_dbl(sim0, ~ .x$N[[1]]),
N_p = map_dbl(sim0, ~ .x$N[[2]])) %>%
select(-sim0, -starts_with("invest")) %>%
pivot_longer(starts_with("N"), names_to = "spp",
values_to = "N") %>%
mutate(spp = ifelse(str_ends(spp, "_c"), "conflicter", "partitioner") %>%
factor(levels = c("conflicter", "partitioner"))) %>%
ggplot(aes(dp, N, color = spp)) +
# geom_vline(xintercept = c(0.3, 0.7), size = 0.5, color = "gray70",
# linetype = "22") +
geom_line(size = 0.5) +
geom_point(size = 1.5) +
scale_color_manual(values = plasma(100)[c(1, 85, 65)][-3]) +
scale_x_continuous(expression(d[p]), breaks = c(1,3,5,7,9)/10) +
scale_y_continuous("Equilibrium abundance", breaks = 2400 + 800 * 0:2,
limits = c(2400, 4305)) +
theme(axis.title = element_blank(),
axis.text = element_blank(),
legend.position = "none",
plot.margin = margin(0,0,0,0)) +
NULL
if (.RESAVE_PLOTS) save_plot(abund_d_p, 1.75, 1,
.name = "fig5_sub-plots/abund_d_p")
#'
#' Simulate invaders into a 2-species community.
#'
one_dx_dp_inv_sim <- function(.dx, .dp, .sim0, .sim_max_t) {
# .dx = 0.6; .dp = 0.6; .sim_max_t = 10e3L
# rm(out, survs, x_vals, O_inv, sim0, raw_sim0, .V0_0, .d_vec,
# .dx, .dp, .sim_max_t)
stopifnot(!inherits(.sim0, "list"))
.d_vec <- c(- .dx, .dp)
O_inv <- with(.sim0, sum(N * exp(- .d_vec[1] * V1^2 - .d_vec[2] * V2^2)))
x_vals <- seq(0.1, 2.5, 0.1) %>%
round(1)
survs <- x_vals %>%
map(function(.x) {
#' No point in doing simulation if we know the invader's fitness
#' will be < 1
f0 <- with(formals(quant_gen), {
exp(r0 - f * .x^2 - a0 - a0 * exp(-.x^2) * O_inv)
})
if (f0 < 1) return(c(rep(1L, nrow(.sim0)), 0L))
.N0 <- c(.sim0$N, 1)
.V0 <- cbind(rbind(.sim0$V1, .sim0$V2), c(.x, 0)) %>%
list()
raw_sim1 <- tibble(eta = 0, d1 = .d_vec[1], d2 = .d_vec[2],
final_t = .sim_max_t, spp_gap_t = 0L,
save_every = 1L,
V0 = .V0, N0 = list(.N0)) %>%
run_determ_1rep()
surv_spp <- raw_sim1 %>%
select(spp, time, N, V1, V2) %>%
unnest(c(spp, time, N, V1, V2)) %>%
filter(time == max(time)) %>%
.[["spp"]] %>%
as.integer()
surv_vec <- as.integer(1:(nrow(.sim0)+1) %in% surv_spp)
# Stable equilibrium?
if (!equil_check(raw_sim1)) {
warning("\nDidn't reach equilibrium at `x` = ", .x,
" for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
surv_vec[surv_vec == 1] <- NA_integer_
}
if ((max(Re(raw_sim1[["eig"]][[1]])) - 1) > 1e-6) {
warning("\nEquilibrium unstable at `x` = ", .x,
" for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
surv_vec[surv_vec == 1] <- NA_integer_
}
return(surv_vec)
})
out <- cbind(x_vals, do.call(rbind, survs)) %>%
as.data.frame() %>%
set_names(c("x", "spp1", "spp2", "spp3")) %>%
as_tibble() %>%
mutate(dx = .dx, dp = .dp) %>%
select(dx, dp, everything())
return(out)
}
# # Takes ~9 min
# dx_dp_sims <- invest_d_sims %>%
# select(dx, dp, sim0) %>%
# rename_with(~ paste0(".", .x)) %>%
# # Some combos need to be run longer to reach equilibrium:
# mutate(.sim_max_t = case_when(.dp >= 0.7 & .dp <= 0.85 & .dx >= 0.7 ~ 30e3L,
# .dp >= 0.7 & .dp <= 0.85 & .dx > 0.6 ~ 20e3L,
# .dp >= 0.7 & .dp <= 0.85 & .dx >= 0.2 ~ 30e3L,
# TRUE ~ 10e3L)) %>%
# # head() %>%
# as.list() %>%
# c(list(FUN = one_dx_dp_inv_sim, SIMPLIFY = FALSE)) %>%
# do.call(what = mcmapply) %>%
# do.call(what = bind_rows)
# saveRDS(dx_dp_sims, "_results/results-rds/dx_dp_sims.rds")
dx_dp_sims <- readRDS("_results/results-rds/dx_dp_sims.rds")
#' This should return no rows bc it checks for any reps not reaching an
#' equilibrium or reach an unstable equilibrium
dx_dp_sims %>%
filter(if_any(starts_with("spp"), is.na))
#' Only 4 possible outcomes:
#' 1. both residents persist, invader excluded
#' 2. species 1 excluded, invader and species 2 coexist
#' 3. all 3 species persist
#' 4. both residents excluded, invader persists
dx_dp_sims %>%
mutate(outcome = paste0(spp1, spp2, spp3)) %>%
distinct(outcome)
dx_dp_sims %>% filter(dx == 0.5 & (dp == 0.3 | dp == 0.7) & x == 1)
invest_d_comm_p <- invest_d_sims %>%
filter(dx == 0.5, dp %in% c(0.3, 0.7)) %>%
arrange(dp) %>%
.[["sim0"]] %>%
map(function(dd) {
dd %>%
ggplot(aes(V1, V2, color = spp))+
geom_abline(slope = 1, intercept = 0,
linetype = 2, color = "gray70") +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_vline(xintercept = 0, color = "gray70", size = 0.5) +
geom_point(shape = 19, size = 3) +
scale_color_manual(values = plasma(100)[c(1, 85)], guide = "none") +
scale_y_continuous("Partitioning investment", breaks = 0.5 * 0:2) +
scale_x_continuous("Conflict investment", breaks = 0.5 * 0:2) +
coord_equal(xlim = c(-0.05, 1.3), ylim = c(-0.05, 1.3),
clip = "off") +
theme(panel.background = element_blank(),
axis.title = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
plot.margin = margin(0,0,0,0)) +
NULL
})
if (.RESAVE_PLOTS) {
save_plot(invest_d_comm_p[[1]], 0.8, 0.8, .name = "fig5_sub-plots/comm_0.3")
save_plot(invest_d_comm_p[[2]], 0.8, 0.8, .name = "fig5_sub-plots/comm_0.7")
}
invade_dp <- function(.dp, .inv_x = 1) {
# .dp = 0.7; .inv_x = 1
# rm(.dx, .dp, .inv_x, equil, raw_sim1, prequel, out)
.dx <- 0.5
stopifnot(length(.dp) == 1 && is.numeric(.dp) && .dp > 0)
stopifnot(length(.inv_x) == 1 && is.numeric(.inv_x) && .inv_x > 0)
equil <- invest_d_sims %>%
filter(dx == .dx, dp == .dp) %>%
.[["sim0"]] %>%
.[[1]]
raw_sim1 <- tibble(eta = 0, d1 = -.dx, d2 = .dp,
final_t = 10e3L, spp_gap_t = 0L,
save_every = 1L,
V0 = list(rbind(c(equil$V1, .inv_x),
c(equil$V2, 0))),
N0 = list(c(equil$N, 1))) %>%
run_determ_1rep() %>%
select(spp, time, N, V1, V2) %>%
unnest(c(spp, time, N, V1, V2)) %>%
mutate(Omega = calc_omega(., d_vec = c(-.dx, .dp)),
invest = sqrt(V1^2 + V2^2),
dp = .dp)
prequel <- equil %>% mutate(dp = .dp, time = -1L)
out <- bind_rows(prequel, raw_sim1) %>%
select(dp, everything())
return(out)
}
invade_dp_df <- map_dfr(c(0.3, 0.7), invade_dp)
invade_dp_p_list <- sort(unique(invade_dp_df$dp)) %>%
map(function(.dp) {
.dd <- invade_dp_df %>% filter(dp == .dp)
ext_df <- .dd %>%
group_by(spp) %>%
filter(time == max(time)) %>%
ungroup() %>%
filter(time < max(time))
y_scales <- list(invest = scale_y_continuous("Total investment",
limits = c(-0.16, 1.6),
breaks = c(0, 0.7, 1.4)),
Omega = scale_y_continuous("ECN",
limits = c(-1100, 11000),
breaks = 5e3*0:2),
N = scale_y_continuous("Abundance",
limits = c(-410, 4100),
breaks = 2e3*0:2))
plots <- imap(y_scales, function(y_scale, y_var) {
.dd %>%
filter(time < 350) %>%
mutate(time = ifelse(time < 0, time * 100, time)) %>%
ggplot(aes(time, !!sym(y_var), color = spp)) +
geom_vline(xintercept = 0, color = "gray70", size = 0.5) +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_line(size = 0.75) +
geom_point(data = ext_df, size = 2, shape = 4) +
scale_color_manual(values = plasma(100)[c(1, 85, 65)],
guide = "none") +
scale_x_continuous("Time", breaks = 100 * 0:3) +
y_scale +
theme(axis.title = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_blank(),
plot.margin = margin(0,0,0,0)) +
coord_cartesian(xlim = c(-50, NA), expand = FALSE)
})
p <- wrap_plots(plots, ncol = 1) &
theme(plot.margin = margin(0,0,0,0))
return(p)
})
# wrap_plots(invade_dp_p_list, ncol = 2)
if (.RESAVE_PLOTS) {
save_plot(invade_dp_p_list[[1]], 2, 1.75, .name = "fig5_sub-plots/invade_0.3")
save_plot(invade_dp_p_list[[2]], 2, 1.75, .name = "fig5_sub-plots/invade_0.7")
}
# =============================================================================*
# =============================================================================*
# Full outcomes
# =============================================================================*
# =============================================================================*
d_outcome_p <- dx_dp_sims %>%
mutate(outcome = case_when(spp1 == 1 & spp3 == 1 ~ "3-spp coexist",
spp1 == 1 & spp3 == 0 ~ "invader excluded",
spp1 == 0 & spp2 == 1 & spp3 == 1 ~ "conflict resident excluded",
spp1 == 0 & spp2 == 0 & spp3 == 1 ~ "both residents excluded",
TRUE ~ NA_character_) %>%
factor(levels = c("both residents excluded",
"conflict resident excluded",
"3-spp coexist",
"invader excluded"))) %>%
mutate(dx_f = factor(dx, levels = sort(unique(dx)),
labels = paste("d[x] ==", sort(unique(dx))))) %>%
ggplot(aes(dp, x)) +
geom_raster(aes(fill = outcome)) +
scale_fill_manual(NULL, values = c("gray80", plasma(100)[c(70, 40, 20)])) +
scale_y_continuous(expression("Invader" ~ x[0]), breaks = 0:2) +
scale_x_continuous(expression(d[p]), breaks = c(2,8,14)/10) +
facet_wrap(~ dx_f, labeller = label_parsed, nrow = 3) +
guides(fill = guide_legend(ncol = 1)) +
theme(legend.position = "right") +
NULL
if (.RESAVE_PLOTS) save_plot(d_outcome_p, 6, 4, .prefix = "S5-", .png = TRUE)
dx_dp_minx <- dx_dp_sims %>%
split(interaction(.$dx, .$dp)) %>%
map_dfr(function(.d) {
.res_excl <- .d %>%
filter(spp1 == 0 & spp3 == 1)
if (nrow(.res_excl) == 0) {
.min_x <- 100
} else {
.min_x <- min(.res_excl$x)
}
out <- select(.d[1,], dx, dp) %>%
mutate(min_x = .min_x)
return(out)
}) %>%
mutate(inf = factor(min_x == 100),
min_x = ifelse(min_x == 100, 2.2, min_x))
dp_minx_p <- dx_dp_minx %>%
filter(dx %in% (c(3, 5, 7, 9) / 10)) %>%
filter(dp <= 0.9) %>%
group_by(dx) %>%
filter(inf == "FALSE" | dp ==
if (any(inf == "TRUE")) max(dp[inf == "TRUE"]) else -Inf) %>%
ungroup() %>%
ggplot(aes(dp, min_x)) +
geom_hline(yintercept = 0.1, color = "gray70", size = 0.5) +
geom_hline(yintercept = 2.2, color = "gray70", size = 0.5) +
geom_line(aes(color = factor(dx)), size = 0.5) +
geom_point(aes(shape = inf, color = factor(dx)), size = 1.5) +
scale_shape_manual(values = c(19, 1), guide = "none") +
# scale_color_viridis_d(expression(d[x]), begin = 0.2, end = 0.8) +
scale_color_manual(values = c("gray80", "gray60", "gray40", "black")) +
scale_x_continuous(expression(d[p]), breaks = c(1,5,9)/10) +
scale_y_continuous(expression(x[inv]), breaks = 0:2) +
theme(axis.title = element_blank(),
axis.text = element_blank(),
legend.position = "none",
plot.margin = margin(0,0,0,0)) +
NULL
if (.RESAVE_PLOTS) save_plot(dp_minx_p, 2, 1.75,
.name = "fig5_sub-plots/dp_minx")
paired_invasions <- function(.dx, .dp, .steps = c(-0.15, 0.15)) {
# .dx = 0.7; .dp = 0.6; .step = 0.15
# rm(.d_vec, .min_x, .x_vec, equil, out_list, out, i, .x, raw_sim1, prequel)
stopifnot(length(.steps) == 2 && is.numeric(.steps) &&
.steps[1] < 0 && .steps[2] > 0)
.d_vec = c(- abs(.dx), abs(.dp))
.min_x <- dx_dp_minx %>%
filter(dx == abs(.d_vec[1]), dp == .d_vec[2]) %>%
.[["min_x"]]
.x_vec <- .min_x + .steps
equil <- invest_d_sims %>%
filter(dx == abs(.d_vec[1]), dp == .d_vec[2]) %>%
.[["sim0"]] %>%
.[[1]]
out_list <- rep(list(NA), 2)
for (i in 1:length(.x_vec)) {
.x = .x_vec[i]
raw_sim1 <- tibble(eta = 0, d1 = .d_vec[1], d2 = .d_vec[2],
final_t = 10e3L, spp_gap_t = 0L,
save_every = 1L,
V0 = list(rbind(c(equil$V1, .x),
c(equil$V2, 0))),
N0 = list(c(equil$N, 1))) %>%
run_determ_1rep()
# Stable equilibrium?
if (!equil_check(raw_sim1)) {
stop("\nDidn't reach equilibrium at `x` = ", .x,
" for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
}
if ((max(Re(raw_sim1[["eig"]][[1]])) - 1) > 1e-6) {
stop("\nEquilibrium unstable at `x` = ", .x,
" for `.d_vec` = ",
paste(.d_vec, collapse = ", "))
}
out <- raw_sim1 %>%
select(spp, time, N, V1, V2) %>%
unnest(c(spp, time, N, V1, V2)) %>%
mutate(Omega = calc_omega(., d_vec = .d_vec),
invest = sqrt(V1^2 + V2^2),
x_inv = .x) %>%
select(x_inv, everything())
prequel <- equil %>% mutate(x_inv = .x, time = -1L)
out_list[[i]] <- bind_rows(prequel, out)
}
do.call(bind_rows, out_list)
}
# The below code shows plots of the weird edge case issue in Figure S5.
# x <- paired_invasions(0.1, 1.35, .steps = c(-2, 0.5)) %>%
# filter(x_inv == min(x_inv))
# y <- paired_invasions(0.3, 1.35, .steps = c(-0.01, 0.4)) %>%
# filter(x_inv == max(x_inv))
# # The problem one:
# z <- paired_invasions(0.2, 1.35, .steps = c(-0.6, 0.9)) %>%
# filter(x_inv == min(x_inv))
#
# problem_plot <- function(.d, .dx, y_var, .max, .leg = TRUE) {
# p <- .d %>%
# ggplot(aes(time, !!sym(y_var), color = spp)) +
# ggtitle(sprintf("dx = %.2f", .dx)) +
# geom_line() +
# scale_y_continuous(limits = c(0, .max)) +
# scale_color_manual(values = plasma(100)[c(1, 85, 65)])
# if (!.leg) p <- p +
# theme(legend.position = "none")
# return(p)
# }
#
#
# problem_plot(x, 0.1, "invest", 1.4) +
# problem_plot(z, 0.2, "invest", 1.4) +
# problem_plot(y, 0.3, "invest", 1.4) +
# problem_plot(x, 0.1, "Omega", 7000) +
# problem_plot(z, 0.2, "Omega", 7000) +
# problem_plot(y, 0.3, "Omega", 7000) +
# problem_plot(x, 0.1, "N", 4700) +
# problem_plot(z, 0.2, "N", 4700) +
# problem_plot(y, 0.3, "N", 4700) +
# plot_layout(ncol = 3, byrow = TRUE, guides = "collect")
paired_inv_df <- paired_invasions(0.5, 0.6, .steps = c(-0.9, 0.1))
paired_inv_p <- sort(unique(paired_inv_df$x_inv)) %>%
map(function(xinv) {
.dd <- paired_inv_df %>% filter(x_inv == xinv)
ext_df <- .dd %>%
group_by(spp) %>%
filter(time == max(time)) %>%
ungroup() %>%
filter(time < max(time))
y_scales <- list(invest = scale_y_continuous("Total investment",
limits = c(-0.17, 1.7),
breaks = c(0, 0.8, 1.6)),
Omega = scale_y_continuous("ECN",
limits = c(-900, 9000),
breaks = 4e3*0:2),
N = scale_y_continuous("Abundance",
limits = c(-410, 4100),
breaks = 2e3*0:2))
plots <- imap(y_scales, function(y_scale, y_var) {
.dd %>%
filter(time < 400) %>%
mutate(time = ifelse(time < 0, time * 100, time)) %>%
ggplot(aes(time, !!sym(y_var), color = spp)) +
geom_vline(xintercept = 0, color = "gray70", size = 0.5) +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_line(size = 0.75) +
geom_point(data = ext_df, size = 2, shape = 4) +
scale_color_manual(values = plasma(100)[c(1, 85, 65)],
guide = "none") +
scale_x_continuous("Time", breaks = 100 * 0:3) +
y_scale +
theme(axis.title = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_blank(),
plot.margin = margin(0,0,0,0)) +
coord_cartesian(xlim = c(-50, NA), expand = FALSE)
})
p <- wrap_plots(plots, ncol = 1) &
theme(plot.margin = margin(0,0,0,0))
return(p)
})
# wrap_plots(paired_inv_p, ncol = 2)
if (.RESAVE_PLOTS) {
save_plot(paired_inv_p[[1]], 2, 1.75, .name = "fig5_sub-plots/under")
save_plot(paired_inv_p[[2]], 2, 1.75, .name = "fig5_sub-plots/over")
}
paired_inv_comm_p <- invest_d_sims %>%
filter(dx == 0.5, dp == 0.6) %>%
.[["sim0"]] %>%
.[[1]] %>%
ggplot(aes(V1, V2, color = spp))+
geom_abline(slope = 1, intercept = 0,
linetype = 2, color = "gray70") +
geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
geom_vline(xintercept = 0, color = "gray70", size = 0.5) +
geom_point(shape = 19, size = 3) +
# geom_point(shape = 19, size = 3, alpha = 0.75) +
# geom_point(shape = 1, size = 3) +
geom_segment(y = -0.05, yend = 0.05,
x = dx_dp_minx %>%
filter(dx == 0.5, dp == 0.6) %>%
.[["min_x"]],
xend = dx_dp_minx %>%
filter(dx == 0.5, dp == 0.6) %>%
.[["min_x"]], color = "black", size = 0.5) +
scale_color_manual(values = plasma(100)[c(1, 85)], guide = "none") +
scale_y_continuous("Partitioning investment", breaks = 0:2) +
scale_x_continuous("Conflict investment", breaks = 0:2) +
coord_equal(xlim = c(-0.025, 2), ylim = c(-0.025, 2), clip = "off") +
theme(panel.background = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
plot.margin = margin(0,0,0,0))
paired_inv_comm_p_invasion <-
list("min", "max") %>%
set_names() %>%
map(~ eval(parse(text = .x))) %>%
map(function(.f) {
.dd <- filter(paired_inv_df, time == 0, spp == 3) %>%
filter(V1 == .f(V1))
p <- paired_inv_comm_p +
geom_point(data = .dd,
shape = 19, size = 3, color = plasma(100)[c(65)])
gt <- ggplot_gtable(ggplot_build(p))
gt$layout$z[grepl("^panel", gt$layout$name)] <- max(gt$layout$z) + 1
gt
})
if (.RESAVE_PLOTS) {
save_plot(paired_inv_comm_p, 1, 1, .name = "fig5_sub-plots/comm")
save_plot(function() {grid.newpage(); grid.draw(paired_inv_comm_p_invasion$max)},
0.8, 0.8, .name = "fig5_sub-plots/comm-above")
save_plot(function() {grid.newpage(); grid.draw(paired_inv_comm_p_invasion$min)},
0.8, 0.8, .name = "fig5_sub-plots/comm-below")
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.