In lucasnell/myvatn_terrestrial_foodweb: Myvatn Terrestrial Foodweb
Appendix S2: Supplementary Figures
knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE,
fig.align = "center")
knitr::opts_knit$set(root.dir = normalizePath(".."))
# load packages
suppressPackageStartupMessages({
library(mtf)
library(tidyverse)
library(grid)
library(parallel)
library(cowplot)
library(egg)
})
# This sets plotting device on LAN computer:
if (!isTRUE(getOption('knitr.in.progress')) && file.exists(".Rprofile")) {
source(".Rprofile")
}
# Clean captions
cc <- function(.x) {
.x <- gsub("\n", "", .x)
.x <- gsub("\\s+", " ", .x)
return(.x)
}
parlist <- par_estimates %>%
filter(V==1, H==1, X==1, iI == 10) %>%
as.list()
V_gain <- function(V, D, aDV = parlist[["aDV"]]) {
hD <- parlist[["hD"]]
(aDV*D*V/(1 + aDV*hD*D)) / V
}
V_loss <- function(V, X, H, M, q, hM = parlist[["hM"]], .no_M_to_X = FALSE) {
aX <- parlist[["aX"]]
hX <- parlist[["hX"]]
if (!.no_M_to_X) {
Vl <- ((aX*V*X)/(1 + aX*hX*(V + H) + (aX * q)*hM*M)) / V
} else {
Vl <- ((aX*V*X)/(1 + aX*hX*(V + H))) / V
}
return(Vl)
}
H_gain <- function(P, H, aPH = parlist[["aPH"]]) {
hP <- parlist[["hP"]]
(aPH*P*H/(1 + aPH*hP*P)) / H
}
H_loss <- function(H, X, V, M, q, hM = parlist[["hM"]], .no_M_to_X = FALSE) {
aX <- parlist[["aX"]]
hX <- parlist[["hX"]]
if (!.no_M_to_X) {
Hl <- ((aX*H*X)/(1 + aX*hX*(V + H) + (aX * q)*hM*M)) / H
} else {
Hl <- (aX*H*X)/(1 + aX*hX*(V + H)) / H
}
return(Hl)
}
upper_levels <- c("detritivore", "herbivore", "predator")
nae_caption <- "Time series of proportional changes in N content among all six
pools with no midge pulse and when starting N content
($N_0$) was randomized.
We ran 20 repetitions, each with different random starting N
for each pool
[$N_0 \\sim \\text{Uniform}(0.5 \\: N_{eq}, \\; 1.5 \\: N_{eq})$].
Sub-panels separate repetitions.
Lines indicate the N content at time $t$ relative to the
equilibrium value [i.e. $(N_t - N_{eq}) / N_{eq}$].
Points indicate the starting N content relative to the
equilibrium value [i.e. $(N_0 - N_{eq}) / N_{eq}$].
Parameter values are as in Table 1."
nae_caption_short <- "Time series for the return to equilibrium N content"
not_at_eq <- function(.r) {
# Equilibrium Ns
Neqs <- par_estimates %>%
filter(V==1, H==1, X==1, iI == 10) %>%
select(ends_with("eq")) %>%
as.list()
rename_vec <- list(I = "soil", D = "detritus", P = "plant", V = "detritivore",
H = "herbivore", X = "predator", M = "midge")
names(Neqs) <- names(Neqs) %>%
gsub(pattern = "eq$", replacement = "") %>%
map_chr(~ rename_vec[[.x]])
# Random starting Ns
N0s <- par_estimates %>%
filter(V==1, H==1, X==1, iI == 10) %>%
select(ends_with("eq")) %>%
mutate_all(~ runif(1, .x * 0.5, .x * 1.5)) %>%
as.list()
names(N0s) <- gsub("eq$", "0", names(N0s))
.dd <- food_web(tmax = 100, s = 1000, b = 0, w = 0, other_pars = list(q = 3),
pool_starts = N0s) %>%
filter(pool != "midge") %>%
mutate(pool = droplevels(pool),
rep = .r) %>%
group_by(pool) %>%
mutate(rel_N = (N - Neqs[[pool[[1]]]]) / Neqs[[pool[[1]]]]) %>%
ungroup()
return(.dd)
}
set.seed(1755669093)
nae_df <- map_dfr(1:20, not_at_eq) %>%
mutate(rep = factor(rep))
nae_p <- nae_df %>%
ggplot(aes(time, rel_N, color = pool)) +
geom_hline(yintercept = 0, color = "gray70") +
geom_point(data = nae_df %>%
filter(time == 0) %>%
mutate(time = -5)) +
geom_line(size = 0.5) +
scale_y_continuous("Proportional change in N", breaks = c(-0.5, 0, 0.5)) +
scale_x_continuous("Time (days)", breaks = c(0, 50, 100)) +
facet_wrap(~ rep, ncol = 4) +
scale_color_manual(NULL, values = color_pal()[c(4:6, 1:3)]) +
theme(panel.spacing = unit(1.5, "lines"),
strip.text = element_blank()) +
NULL
nae_p
\clearpage
supp_sims1 <- crossing(.aDV = c(par_estimates$aDV[1], par_estimates$aPH[1]),
.aPH = c(par_estimates$aPH[1], par_estimates$aDV[1])) %>%
pmap_dfr(function(.aDV, .aPH) {
food_web(tmax = 100, s = 10, b = 20, w = 20,
other_pars = list(q = 3, aDV = .aDV, aPH = .aPH)) %>%
mutate(aDV = .aDV, aPH = .aPH)
}) %>%
spread(pool, N) %>%
mutate(Vg = V_gain(detritivore, detritus, aDV = aDV),
Hg = H_gain(plant, herbivore, aPH = aPH)) %>%
select(-soil:-midge) %>%
gather("variable", "value", Vg:Hg) %>%
mutate(pool = factor(ifelse(grepl("^V", variable), "detritivore", "herbivore")),
aDV = factor(aDV, levels = range(as.numeric(paste(aDV))),
labels = paste(c("low", "high"), "detritivore uptake")),
aPH = factor(aPH, levels = range(as.numeric(paste(aPH))),
labels = paste(c("low", "high"), "herbivore uptake"))) %>%
select(-variable) %>%
select(aDV, aPH, pool, everything()) %>%
group_by(aDV, aPH, pool) %>%
mutate(value = value - value[1]) %>%
ungroup()
s1_caption <- "Time series of the bottom-up
effects on detritivore and herbivore pools with (a,b) low and
(c,d) high herbivore uptake rates ($a_P$), and with (a,c) low and
(b,d) high detritivore uptake rates ($a_D$).
Parameter values are as in Table 1, with low $a_D = 0.012$,
low $a_P = 0.012$, high $a_D = 0.032$, high $a_P = 0.032$,
$w = 20$, $b = 20$, and $q = 3$."
s1_caption_short <- "Effects of varying herbivore and detritivore uptake rates on
bottom-up effects"
trans_p3 <- ggplot(data = NULL) +
geom_hline(yintercept = 0, color = "gray70") +
geom_line(data = supp_sims1,
aes(time, value,
color = pool), size = 1, linetype = 1) +
scale_y_continuous(expression("Bottom-up effect on pool (" * day^{-1} * ")" )) +
scale_x_continuous("Time (days)") +
scale_color_manual(values = color_pal()[1:2]) +
geom_text(data = tibble(time = c( 35, 28),
value = c(0.070, 0.007),
aDV = sort(unique(supp_sims1$aDV))[1],
aPH = sort(unique(supp_sims1$aPH))[1],
lab = c("detritivore", "herbivore")),
aes(time, value, label = lab), color = color_pal()[1:2],
hjust = 0, vjust = 0, lineheight = 0.75, size = 10 / 2.835, parse = TRUE) +
facet_grid(aPH ~ aDV) +
theme(legend.position = "none",
strip.text = element_text(face = "plain", size = 11,
margin = margin(b = 4, l = 4)),
panel.spacing.x = unit(2.5, "lines"),
strip.background = element_blank(),
axis.text = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12)) +
NULL
trans_p3
vain_equil_caption <- "Equilibrium N pool size for all six
pools versus the inorganic nutrient input rate ($i_I$).
Line type indicates the density-independent loss of soil
nutrients ($\\mu_I$).
Parameter values are as in Table 1."
vain_equil_caption_short <- "Effects of the inorganic nutrient input rate and
density-independent loss of soil
nutrients on equilibrium N pool sizes"
vary_iI_muI_equil <- function(.) {
Z <- equil_pools(iI = .$iI[1], muI = .$muI[1])
mutate(Z[["vals"]], iI = .$iI[1], muI = .$muI[1])
}
# Takes ~12 sec with 4 threads
vain_equil_df <- crossing(iI = seq(0, 20, length.out = 61),
muI = c(0.005, 0.01, 0.02)) %>%
split(1:nrow(.)) %>%
mclapply(vary_iI_muI_equil,
mc.cores = parallel::detectCores()) %>%
bind_rows()
# This just sets the bounds for the facets:
vain_dummy_equil <- tibble(iI = 20,
pool = factor(c("soil", "detritus", "plant", upper_levels),
levels = c("soil", "detritus", "plant", upper_levels)),
eq = c(700, 60, 60, 2, 2, 2))
vain_equil_p <- vain_equil_df %>%
mutate_at(vars(muI), factor) %>%
ggplot(aes(iI, eq)) +
geom_hline(yintercept = 0, color = "gray70") +
geom_line(aes(color = pool, linetype = muI), size = 0.75) +
geom_blank(data = vain_dummy_equil) +
scale_y_continuous(expression("Equilibrium pool size (g N" ~ m^{-2} * ")")) +
scale_x_continuous(expression(paste("Inorganic nutrient input rate",
~ "(" * g ~ N ~ m^{-2} ~ d^{-1} * ")"))) +
facet_wrap(~ pool, nrow = 2, scales = "free_y") +
scale_color_manual(NULL, values = color_pal()[c(4:6, 1:3)], guide = FALSE) +
scale_linetype_manual("Density-independent\nloss of soil nutrients:", values = c(2, 1, 3)) +
theme(panel.spacing = unit(1, "lines"), legend.position = "top") +
NULL
vain_equil_p
vain_ts_caption <- "Time series of proportional changes in N content among all six
pools with varying values of $i_I$ and $\\mu_I$,
the inorganic nutrient input rate and density-independent
loss of soil nutrients, respectively.
Sub-panel labels indicate the values of each parameter, and
middle values are the defaults used in the main text.
Lines indicate the N content at time $t$ relative to the
starting value [i.e. $(N_t - N_{0}) / N_{0}$].
Parameter values are as in Table 1, with
predator exploitation $q = 3$, pulse duration $w = 0$, and
pulse rate $b = 0$.
Note that when $i_I = 5$ and $\\mu_I = 0.02$, predators
go extinct so are not shown in that panel."
vain_ts_caption_short <- "Time series of N content in response to a midge pulse,
while varying the inorganic nutrient input rate and
density-independent loss of soil nutrients"
vary_iI_muI <- function(.iI, .muI) {
# .iI <- 5; .muI <- 0.02
# rm(.iI, .muI, pools, sim)
if (.iI == 5 && .muI == 0.02) {
pools <- equil_pools(iI = .iI, muI = .muI, Xeq = 0)
sim <- food_web(tmax = 100, s = 10, b = 50, w = 20,
other_pars = list(q = 3, iI = .iI, muI = .muI),
ep_obj = pools) %>%
filter(pool != "midge") %>%
mutate(pool = droplevels(pool)) %>%
filter(pool != "predator")
} else {
pools <- equil_pools(iI = .iI, muI = .muI)
sim <- food_web(tmax = 100, s = 10, b = 50, w = 20,
other_pars = list(q = 3, iI = .iI, muI = .muI),
ep_obj = pools) %>%
filter(pool != "midge") %>%
mutate(pool = droplevels(pool))
}
sim <- sim %>%
group_by(pool) %>%
mutate(N_rel = (N - N[1]) / (N[1])) %>%
ungroup() %>%
mutate(iI = .iI, muI = .muI)
return(sim)
}
vain_df <- crossing(.iI = c(5, 10, 20),
.muI = c(0.005, 0.01, 0.02)) %>%
pmap_dfr(vary_iI_muI) %>%
mutate_at(vars(iI, muI), factor) %>%
mutate(muI = factor(muI, levels = levels(muI),
labels = paste("italic(mu[I]) ==", levels(muI))),
iI = factor(iI, levels = levels(iI),
labels = paste("italic(i[I]) ==", levels(iI))))
vain_dummy <- vain_df %>%
group_by(iI, muI) %>%
filter(N_rel == max(N_rel)) %>%
ungroup() %>%
select(-pool, -N) %>%
mutate(time = 10,
time2 = 10 + 20,
N_rel = -0.1 * N_rel)
vain_p <- vain_df %>%
ggplot(aes(time, N_rel)) +
geom_hline(yintercept = 0, color = "gray70") +
geom_segment(data = vain_dummy,
aes(xend = time2, yend = N_rel), size = 1.5) +
geom_line(aes(color = pool), size = 0.5) +
scale_y_continuous("Proportional change in N") +
scale_x_continuous("Time (days)") +
facet_wrap(iI ~ muI, nrow = 3, scales = "free_y", labeller = label_parsed) +
scale_color_manual(NULL, values = color_pal()[c(4:6, 1:3)]) +
guides(color = guide_legend(override.aes = list(size = 1))) +
theme(panel.spacing = unit(1, "lines"),
strip.text = element_text(family = "CMU Serif")) +
NULL
vain_p
fig2_abs_caption <- "Time series of N content among all six pools.
Lines indicate the N content at time $t$.
The lower black bar represents the duration of the midge N pulse.
Parameter values are as in Table 1, with
predator exploitation $q = 3$,
pulse duration $w = 20$, and pulse rate $b = 50$."
fig2_abs_caption_short <- "Time series of absolute N content in response to a midge pulse"
middle_sim <- food_web(tmax = 100, s = 10, b = 50, w = 20,
other_pars = list(q = 3))
# This just sets the bounds for the facets:
dummy_data <- tibble(time = 20,
pool = factor(rep(c("soil", "detritus", "plant", upper_levels),
each = 2),
levels = c("soil", "detritus", "plant", upper_levels)),
N = c(middle_sim %>%
filter(pool == "soil") %>%
.[["N"]] %>%
max() %>%
rep(2),
middle_sim %>%
filter(pool %in% c("detritus", "plant")) %>%
.[["N"]] %>%
max() %>%
rep(2 * 2),
middle_sim %>%
filter(pool %in% upper_levels) %>%
.[["N"]] %>%
max() %>%
rep(3 * 2))) %>%
group_by(pool) %>%
mutate(N = N[1] * c(-0.1, 1)) %>%
ungroup()
trans_pS1 <- middle_sim %>%
filter(pool != "midge") %>%
mutate(pool = droplevels(pool)) %>%
ggplot(aes(time, N)) +
geom_hline(yintercept = 0, color = "gray70") +
geom_segment(data = tibble(time = 10, time2 = 10+20, N = 0),
aes(xend = time2, yend = N), size = 1.5) +
geom_line(aes(color = pool), size = 1) +
geom_blank(data = dummy_data) +
scale_y_continuous(expression("Nitrogen" ~ (g ~ N ~ m^{-2}))) +
scale_x_continuous("Time (days)") +
facet_wrap(~ pool, nrow = 2, scales = "free") +
scale_color_manual(values = color_pal()[c(4:6, 1:3)], guide = FALSE) +
theme(panel.spacing = unit(1.5, "lines")) +
geom_text(data = dummy_data %>%
filter(N < 0, pool == "soil"),
label = "pulse", size = 10 / 2.835, hjust = 0.5, vjust = 0,
color = "black") +
NULL
trans_pS1
\clearpage
nohv_fig3_caption <- "Time series of proportional changes in N content for
detritivore and herbivore pools
(indicated by the left y-axis).
This is shown for when (a) no detritivores are present,
(b) both are present, and
(c) no herbivores are present.
(a,c) When either detritivores or herbivores were absent,
we solved for the other equilibrium pool sizes
before running the simulations.
The blue shaded regions represent the proportional change in
N for predators, which is indicated by the right y-axis.
Different y-axis scales are used to aid visualization of
the transient dynamics of the herbivores and detritivores,
which had a weaker response to the pulse than predators for
the selected parameter values.
The lower black bar represents the duration of the
midge N pulse.
Parameter values are as in Table 1, with
predator exploitation $q = 3$,
pulse duration $w = 20$, and pulse rate $b = 20$."
nohv_fig3_caption_short <- "Time series of N content in response to a midge
pulse, with herbivores or detritivores removed"
do_nohv_fig3 <- function(.H, .V) {
stopifnot(length(.H) == 1 && .H %in% 0L:1L)
stopifnot(length(.V) == 1 && .V %in% 0L:1L)
if (.H == 0 && .V == 0) {
pools <- equil_pools(Heq = 0, Veq = 0)
}
if (.H == 0 && .V == 1) {
pools <- equil_pools(Heq = 0)
}
if (.H == 1 && .V == 0) {
pools <- equil_pools(Veq = 0)
}
if (.H == 1 && .V == 1) {
pools <- NULL
}
q_ <- 3
sim <- food_web(tmax = 100, s = 10, b = 20, w = 20,
other_pars = list(q = q_),
ep_obj = pools) %>%
filter(pool %in% c(upper_levels, "midge")) %>%
mutate(pool = droplevels(pool)) %>%
mutate(pool = factor(paste(pool), levels = c(upper_levels, "midge")),
H = .H, V = .V)
return(sim)
}
nohv_fig3_df <- crossing(.H = 0L:1L, .V = 0L:1L) %>%
filter(!(.H == 0 & .V == 0)) %>%
pmap_dfr(do_nohv_fig3) %>%
group_by(H, V, pool) %>%
mutate(N_rel = (N - N[1]) / N[1]) %>%
mutate(N_rel = ifelse(is.nan(N_rel), 0, N_rel)) %>%
ungroup()
nohv_fig3_plot <- function(.H, .V,
.mult = 25.2,
.ylims = c(-0.3011, 1.1)) {
# .H = 1; .V = 0
# .mult = 25.2; .ylims = c(-0.3011, 1.1)
# rm(.H, .V, .mult, .ylims)
.title <- case_when(.H == 0 & .V == 1 ~ "No herbivores",
.H == 1 & .V == 0 ~ "No detritivores",
.H == 1 & .V == 1 ~ "Both present",
TRUE ~ "ERROR")
hv_combo <- nohv_fig3_df %>%
filter(H == .H, V == .V) %>%
select(-H, -V)
hv_plot <- hv_combo %>%
filter(!pool %in% c("midge", "predator")) %>%
mutate(pool = droplevels(pool)) %>%
ggplot(aes(time, N_rel)) +
geom_hline(yintercept = 0, color = "gray70") +
geom_segment(data = tibble(time = 10, time2 = 10+20, N_rel = -0.1),
aes(xend = time2, yend = N_rel), size = 1.5) +
geom_ribbon(data = hv_combo %>%
filter(pool == "predator") %>%
mutate(N_rel = N_rel / .mult),
aes(ymin = 0, ymax = N_rel), fill = color_pal(0.5)[3]) +
geom_line(aes(color = pool), size = 0.75) +
scale_x_continuous("Time (days)") +
scale_color_manual(NULL, values = color_pal()) +
ggtitle(.title) +
theme(legend.position = "none",
strip.text.y = element_text(face = "plain", size = 10, angle = 270,
margin = margin(l = 4)),
strip.text.x = element_text(face = "plain", size = 10,
margin = margin(b = 4)),
panel.spacing.y = unit(0, "lines"),
panel.spacing.x = unit(1.5, "lines"),
plot.margin = margin(0,0,t=12,b=12),
plot.title = element_text(size = 12, hjust = 0.5, face = "plain",
margin = margin(0,0,0,b = 10))) +
scale_y_continuous("Proportional change in N", limits = .ylims,
breaks = c(0, 0.5, 1),
sec.axis = sec_axis(~ . * .mult,
name = paste("Proportional change in N",
"(predator)"),
breaks = c(0, 10, 20))) +
NULL
if (.H == 1 && .V == 1) {
hv_plot <- hv_plot +
geom_text(data = tibble(pool = factor(upper_levels[upper_levels !=
"predator"]),
time = c(92, 100),
N_rel = c(0.3, -0.1),
q = factor(paste(rep("high", 2), "exploitation"),
levels = paste(c("low", "high"),
"exploitation"))),
aes(label = pool, color = pool), hjust = 1, size = 10 / 2.835) +
geom_text(data = tibble(time = 20, N_rel = -0.15,
q = factor("high exploitation",
levels = paste(c("low", "high"),
"exploitation"))),
label = "pulse", size = 10 / 2.835, hjust = 0.5, vjust = 1,
color = "black") +
theme(axis.title.y.left = element_text(margin = margin(0,0,0,r=12),
size = 11),
axis.title.y.right = element_text(margin = margin(0,0,0,l=12),
size = 11))
} else {
hv_plot <- hv_plot +
theme(axis.title.y = element_blank())
}
if (.H == 1 && .V == 0) {
hv_plot <- hv_plot +
geom_text(data = tibble(time = 65,
N_rel = 0.8,
q = factor("high exploitation",
levels = paste(c("low", "high"),
"exploitation"))),
label = "predator", color = color_pal()[3], hjust = 1, vjust = 1,
size = 10 / 2.835)
}
if (!(.H == 0 & .V == 1)) {
hv_plot <- hv_plot +
theme(axis.title.x = element_blank(), axis.text.x = element_blank())
} else {
hv_plot <- hv_plot +
theme(axis.title.x = element_text(margin = margin(0,0,0,t=4)))
}
return(hv_plot)
}
nohv_fig3_plot_list <- crossing(.H = 0L:1L, .V = 0L:1L) %>%
filter(!(.H == 0 & .V == 0)) %>%
.[c(2, 3, 1),] %>%
pmap(nohv_fig3_plot) %>%
map(~ .x + theme(axis.line.y.left = element_line(size = 1)))
ggarrange(plots = nohv_fig3_plot_list, ncol = 1, labels = sprintf("(%s)", letters[1:3]),
label.args = list(gp = gpar(fontsize = 14, fontface = "plain"),
vjust = 2, hjust = -1.5))
nohv_fig4_caption <- "Time series of the top-down (\"TD\") and bottom-up (\"BU\")
effects on detritivore and herbivore pools.
This is shown for when (a) no detritivores are present,
(b) both are present, and
(c) no herbivores are present.
(a,c) When either detritivores or herbivores were absent,
we solved for the other equilibrium pool sizes
before running the simulations.
The gray line indicates the top-down effects on both
detritivore and herbivore pools, as they are identical
in the model.
Parameter values are as in Table 1, with
predator exploitation $q = 3$,
pulse duration $w = 20$, and pulse rate $b = 20$."
nohv_fig4_caption_short <- "Time series of top-down and bottom-up effects, with
herbivores or detritivores removed"
do_nohv_fig4 <- function(.H, .V) {
stopifnot(length(.H) == 1 && .H %in% 0L:1L)
stopifnot(length(.V) == 1 && .V %in% 0L:1L)
if (.H == 0 && .V == 0) {
pools <- equil_pools(Heq = 0, Veq = 0)
}
if (.H == 0 && .V == 1) {
pools <- equil_pools(Heq = 0)
}
if (.H == 1 && .V == 0) {
pools <- equil_pools(Veq = 0)
}
if (.H == 1 && .V == 1) {
pools <- NULL
}
q_ <- 3
sim <- food_web(tmax = 100, s = 10, b = 20, w = 20,
other_pars = list(q = q_),
ep_obj = pools) %>%
spread(pool, N) %>%
mutate(Vg = V_gain(detritivore, detritus),
Vl = V_loss(detritivore, predator, herbivore, midge, q_),
Hg = H_gain(plant, herbivore),
Hl = H_loss(herbivore, predator, detritivore, midge, q_)) %>%
select(-soil:-midge) %>%
gather("variable", "value", Vg:Hl) %>%
mutate(pool = factor(ifelse(grepl("^V", variable), "detritivore", "herbivore")),
type = factor(ifelse(grepl("l$", variable), "top-down", "bottom-up"))) %>%
select(-variable) %>%
select(pool, type, everything()) %>%
group_by(pool, type) %>%
mutate(value = value - value[1]) %>%
ungroup() %>%
mutate(H = .H, V = .V)
return(sim)
}
nohv_fig4_df <- crossing(.H = 0L:1L, .V = 0L:1L) %>%
filter(!(.H == 0 & .V == 0)) %>%
pmap_dfr(do_nohv_fig4)
nohv_fig4_plot <- function(.H, .V, .ylims = c(-0.037, 0.144)) {
.title <- case_when(.H == 0 & .V == 1 ~ "No herbivores",
.H == 1 & .V == 0 ~ "No detritivores",
.H == 1 & .V == 1 ~ "Both present",
TRUE ~ "ERROR")
.td_pool <- ifelse(.V == 0, "herbivore", "detritivore")
hv_combo <- nohv_fig4_df %>%
filter(H == .H, V == .V) %>%
select(-H, -V) %>%
mutate(value = ifelse(is.nan(value), 0, value))
hv_plot <- ggplot(data = NULL) +
geom_hline(yintercept = 0, color = "gray70") +
geom_line(data = hv_combo %>% filter(type == "bottom-up"),
aes(time, value,
color = pool,
group = interaction(pool, type)), size = 1, linetype = 1) +
geom_line(data = hv_combo %>% filter(type == "top-down", pool == .td_pool),
aes(time, value), size = 1, color = "gray60") +
xlab("Time (days)") +
ggtitle(.title) +
scale_color_manual(values = color_pal()[1:2]) +
scale_y_continuous(expression("Effect on pool (" * day^{-1} * ")" ),
breaks = c(0, 0.06, 0.12),
limits = .ylims) +
theme(legend.position = "none",
strip.text = element_text(face = "plain", size = 11,
margin = margin(b = 4)),
panel.spacing.x = unit(1.5, "lines"),
plot.margin = margin(0,t=12,b=12,r=4),
strip.background = element_blank(),
axis.text = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12),
plot.title = element_text(size = 12, hjust = 0.5, face = "plain",
margin = margin(0,0,0,b=6)))
if (.H == 1 && .V == 1) {
hv_plot <- hv_plot +
geom_text(data = tibble(time = c( 50),
value = c(0.05),
lab = sprintf("italic(%s)",
c("'BU'['detritivore']"))),
aes(time, value, label = lab),
color = c(color_pal()[1]),
hjust = 0, vjust = 0, lineheight = 0.75,
size = 10 / 2.835, parse = TRUE) +
theme(axis.title.y = element_text(margin = margin(0,0,r=4,l=10)))
} else {
hv_plot <- hv_plot +
theme(axis.title.y = element_blank())
}
if (.H == 1 && .V == 0) {
hv_plot <- hv_plot +
geom_text(data = tibble(time = c(5, 50),
value = c(0.034, 0.05),
lab = sprintf("italic(%s)",
c("'BU'['herbivore']",
"'TD'['both']"))),
aes(time, value, label = lab),
color = c(color_pal()[2], "gray60"),
hjust = 0, vjust = 0, lineheight = 0.75,
size = 10 / 2.835, parse = TRUE)
}
if (!(.H == 0 & .V == 1)) {
hv_plot <- hv_plot +
theme(axis.title.x = element_blank(), axis.text.x = element_blank())
} else {
hv_plot <- hv_plot +
theme(axis.title.x = element_text(margin = margin(0,0,0,t=4)))
}
return(hv_plot)
}
nohv_fig4_plot_list <- crossing(.H = 0L:1L, .V = 0L:1L) %>%
filter(!(.H == 0 & .V == 0)) %>%
.[c(2, 3, 1),] %>%
pmap(nohv_fig4_plot) %>%
map(~ .x + theme(axis.line.y.left = element_line(size = 1)))
ggarrange(plots = nohv_fig4_plot_list, ncol = 1, labels = sprintf("(%s)", letters[1:3]),
label.args = list(gp = gpar(fontsize = 14, fontface = "plain"),
vjust = 2, hjust = -2))
mtop_fig5_caption <- "Cumulative midge effect on
total bottom-up control on detritivores (${}^V BU_{total}$)
total bottom-up control on herbivores (${}^H BU_{total}$),
net top-down control ($TD_{net}$),
top-down intensification, and top-down alleviation
as a function of total midge inputs.
This is shown for when (a) midges only go to the detritus pool,
(b) midges go to both detritus and predator pools, and
(c) midges only go to the predator pool.
The panels show different combinations of high and
low predator exploitation and handling times on midges.
Parameter values are as in Table 1, with
low exploitation $q = 0.8$,
high exploitation $q = 8$,
low midge handling time $h_M = 1.3$,
high midge handling time $h_M = 5.3$,
pulse duration $w = 20$, and
pulse rate $b$ ranging from 0.1 to 23 to produce
different total midge inputs."
mtop_fig5_caption_short <- "Total effect of midges on bottom-up and top-down control,
with varying midge exploitation and handling time,
and with midges going only to detritus or to predators"
do_mtop_fig5 <- function(.midges_not_to) {
# .midges_not_to = "none"
# rm(.midges_not_to)
stopifnot(length(.midges_not_to) == 1 && .midges_not_to %in% c("X", "none", "D"))
.one_combo <- function(.q, .hM, .b) {
# .q = 0.8; .hM = par_estimates$hM[1] * 0.5; .b = 0.1
# rm(.q, .hM, .b)
fw <- food_web(tmax = 500, s = 10, b = .b, w = 20,
other_pars = list(hM = .hM, q = .q),
midges_not_to = .midges_not_to)
fw <- fw %>%
spread(pool, N) %>%
mutate(Vg = V_gain(detritivore, detritus),
Vl = V_loss(detritivore, predator, herbivore, midge, q = .q, hM = .hM,
.no_M_to_X = (.midges_not_to == "X")),
Hg = H_gain(plant, herbivore),
Hl = H_loss(herbivore, predator, detritivore, midge, q = .q, hM = .hM,
.no_M_to_X = (.midges_not_to == "X"))) %>%
summarize(cum_pos_loss_V = sum(Vl[Vl > Vl[1]] - Vl[1]),
cum_pos_loss_H = sum(Hl[Hl > Hl[1]] - Hl[1]),
cum_neg_loss_V = sum(Vl[Vl < Vl[1]] - Vl[1]),
cum_neg_loss_H = sum(Hl[Hl < Hl[1]] - Hl[1]),
cum_gain_V = sum(Vg[Vg > Vg[1]] - Vg[1]),
cum_gain_H = sum(Hg[Hg > Hg[1]] - Hg[1])) %>%
mutate(area = 20 * .b,
q = .q,
hM = .hM) %>%
select(area, q, hM, everything()) %>%
mutate(cum_loss_V = cum_pos_loss_V + cum_neg_loss_V,
cum_loss_H = cum_pos_loss_H + cum_neg_loss_H,
# Changing to magnitudes
cum_neg_loss_V = abs(cum_neg_loss_V),
cum_neg_loss_H = abs(cum_neg_loss_H)) %>%
gather("variable", "value", cum_pos_loss_V:cum_loss_H) %>%
mutate(pool = factor(ifelse(grepl("V$", variable),
"detritivore", "herbivore")),
direction = factor(ifelse(grepl("loss", variable),
"top-down", "bottom-up")),
type = case_when(grepl("^cum_pos_loss", variable) ~
"TD intensification",
grepl("^cum_neg_loss", variable) ~
"TD alleviation",
grepl("^cum_loss_", variable) ~ "TD net",
grepl("^cum_gain_", variable) ~ "BU total",
TRUE ~ NA_character_) %>%
factor(levels = c("BU total", "TD net",
"TD intensification", "TD alleviation"))) %>%
select(-variable)
if (sum(is.na(fw$type)) > 0) stop("Error in processing `type` column.")
return(fw)
}
out <- crossing(.q = c(0.8, 8),
.hM = par_estimates$hM[1] * c(0.5, 2),
.b = seq(0.1, 23, length.out = 21)) %>%
pmap_dfr(.one_combo) %>%
mutate(q = factor(q, levels = sort(unique(q)),
labels = sprintf("%s midge\nexploitation",
c("low", "high"))),
hM = factor(hM, levels = sort(unique(hM)),
labels = paste(c("low", "high"),
"midge\nhandling time"))) %>%
mutate(not_to = .midges_not_to)
return(out)
}
# Takes ~24 sec using 3 threads
mtop_fig5_df <- mclapply(c("X", "none", "D"), do_mtop_fig5,
mc.cores = min(3, parallel::detectCores())) %>%
bind_rows()
mtop_fig5_plot <- function(.midges_not_to,
.ylims = c(-0.688, 4.305)) {
fig5_labs <- list(bquote(italic({}^'V' * 'BU'['total'])),
bquote(italic({}^'H' * 'BU'['total'])),
bquote(italic('TD'['net'])),
bquote(italic('TD'['intensification'])),
bquote(italic('TD'['alleviation'])))
fig5_lvls <- c("detritivore BU total",
"herbivore BU total",
"TD net",
"TD intensification",
"TD alleviation")
.colors <- c(color_pal()[1:2], "gray60", "gray30", "gray60")
# .midges_not_to = "D"; .ylims = c(-0.688, 4.305)
# rm(.midges_not_to, .ylims, .pool)
.title <- case_when(.midges_not_to == "none" ~ "Midges to detritus and predators",
.midges_not_to == "D" ~ "Midges to predators only",
.midges_not_to == "X" ~ "Midges to detritus only",
TRUE ~ "ERROR")
mt_combo <- mtop_fig5_df %>%
filter(not_to == .midges_not_to) %>%
select(-not_to) %>%
mutate(id = ifelse(grepl("^TD", type), paste(type),
sprintf("%s %s", paste(pool), paste(type))),
id = factor(id, levels = fig5_lvls)) %>%
group_by(id, area, q, hM) %>%
summarize(value = value[1]) %>%
ungroup()
mt_plot <- mt_combo %>%
ggplot(aes(area, value)) +
geom_hline(yintercept = 0, color = "gray80") +
geom_line(aes(color = id, linetype = id), size = 0.75) +
ggtitle(.title) +
facet_grid(hM ~ q) +
scale_color_manual(NULL, values = .colors, guide = FALSE) +
scale_linetype_manual(NULL, values = c(1, 1, 1, 3:2),
breaks = fig5_lvls,
labels = fig5_labs) +
guides(linetype = guide_legend(keywidth = 2, ncol = 1, keyheight = 0.6,
override.aes = list(color = .colors, size = 1),
byrow = TRUE)) +
xlab(expression("Total midge input (" * g ~ N ~ m^{-2} * ")")) +
scale_y_continuous("Total midge effect", limits = .ylims, breaks = c(0, 2, 4)) +
theme(legend.background = element_blank(),
legend.position = "right",
legend.direction = "vertical",
legend.justification = "left",
legend.box = "vertical",
legend.title = element_text(size = 11, margin = margin(0,0,0,b=-4)),
legend.text = element_text(size = 10, family = "CMU Serif"),
strip.background = element_blank(),
strip.text.y = element_text(vjust = 0.5, hjust = 0, size = 9, angle = 0),
strip.text.x = element_text(size = 9),
panel.spacing.x = unit(1.5, "lines"),
plot.margin = margin(0,t=12,b=12,r=4),
plot.title = element_text(size = 12, hjust = 0.5, face = "plain",
margin = margin(0,0,0,b=6))) +
NULL
if (.midges_not_to == "none") {
mt_plot <- mt_plot +
theme(axis.title.y = element_text(margin = margin(0,0,0,r=12), size = 12))
} else {
mt_plot <- mt_plot +
theme(axis.title.y = element_blank(),
legend.position = "none")
}
if (.midges_not_to != "D") {
mt_plot <- mt_plot +
theme(axis.title.x = element_blank(), axis.text.x = element_blank())
} else {
mt_plot <- mt_plot +
theme(axis.title.x = element_text(margin = margin(0,0,0,t=4), size = 12))
}
return(mt_plot)
}
mtop_fig5_plot_list <- tibble(.midges_not_to = c("X", "none", "D")) %>%
pmap(mtop_fig5_plot)
ggarrange(plots = mtop_fig5_plot_list, ncol = 1, labels = sprintf("(%s)", letters[1:3]),
label.args = list(gp = gpar(fontsize = 14, fontface = "plain"),
vjust = 2, hjust = -1))
lucasnell/myvatn_terrestrial_foodweb documentation built on Aug. 12, 2020, 6:16 p.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.
Appendix S2: Supplementary Figures
knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE, fig.align = "center") knitr::opts_knit$set(root.dir = normalizePath("..")) # load packages suppressPackageStartupMessages({ library(mtf) library(tidyverse) library(grid) library(parallel) library(cowplot) library(egg) }) # This sets plotting device on LAN computer: if (!isTRUE(getOption('knitr.in.progress')) && file.exists(".Rprofile")) { source(".Rprofile") } # Clean captions cc <- function(.x) { .x <- gsub("\n", "", .x) .x <- gsub("\\s+", " ", .x) return(.x) } parlist <- par_estimates %>% filter(V==1, H==1, X==1, iI == 10) %>% as.list() V_gain <- function(V, D, aDV = parlist[["aDV"]]) { hD <- parlist[["hD"]] (aDV*D*V/(1 + aDV*hD*D)) / V } V_loss <- function(V, X, H, M, q, hM = parlist[["hM"]], .no_M_to_X = FALSE) { aX <- parlist[["aX"]] hX <- parlist[["hX"]] if (!.no_M_to_X) { Vl <- ((aX*V*X)/(1 + aX*hX*(V + H) + (aX * q)*hM*M)) / V } else { Vl <- ((aX*V*X)/(1 + aX*hX*(V + H))) / V } return(Vl) } H_gain <- function(P, H, aPH = parlist[["aPH"]]) { hP <- parlist[["hP"]] (aPH*P*H/(1 + aPH*hP*P)) / H } H_loss <- function(H, X, V, M, q, hM = parlist[["hM"]], .no_M_to_X = FALSE) { aX <- parlist[["aX"]] hX <- parlist[["hX"]] if (!.no_M_to_X) { Hl <- ((aX*H*X)/(1 + aX*hX*(V + H) + (aX * q)*hM*M)) / H } else { Hl <- (aX*H*X)/(1 + aX*hX*(V + H)) / H } return(Hl) } upper_levels <- c("detritivore", "herbivore", "predator")
nae_caption <- "Time series of proportional changes in N content among all six pools with no midge pulse and when starting N content ($N_0$) was randomized. We ran 20 repetitions, each with different random starting N for each pool [$N_0 \\sim \\text{Uniform}(0.5 \\: N_{eq}, \\; 1.5 \\: N_{eq})$]. Sub-panels separate repetitions. Lines indicate the N content at time $t$ relative to the equilibrium value [i.e. $(N_t - N_{eq}) / N_{eq}$]. Points indicate the starting N content relative to the equilibrium value [i.e. $(N_0 - N_{eq}) / N_{eq}$]. Parameter values are as in Table 1." nae_caption_short <- "Time series for the return to equilibrium N content"
not_at_eq <- function(.r) { # Equilibrium Ns Neqs <- par_estimates %>% filter(V==1, H==1, X==1, iI == 10) %>% select(ends_with("eq")) %>% as.list() rename_vec <- list(I = "soil", D = "detritus", P = "plant", V = "detritivore", H = "herbivore", X = "predator", M = "midge") names(Neqs) <- names(Neqs) %>% gsub(pattern = "eq$", replacement = "") %>% map_chr(~ rename_vec[[.x]]) # Random starting Ns N0s <- par_estimates %>% filter(V==1, H==1, X==1, iI == 10) %>% select(ends_with("eq")) %>% mutate_all(~ runif(1, .x * 0.5, .x * 1.5)) %>% as.list() names(N0s) <- gsub("eq$", "0", names(N0s)) .dd <- food_web(tmax = 100, s = 1000, b = 0, w = 0, other_pars = list(q = 3), pool_starts = N0s) %>% filter(pool != "midge") %>% mutate(pool = droplevels(pool), rep = .r) %>% group_by(pool) %>% mutate(rel_N = (N - Neqs[[pool[[1]]]]) / Neqs[[pool[[1]]]]) %>% ungroup() return(.dd) } set.seed(1755669093) nae_df <- map_dfr(1:20, not_at_eq) %>% mutate(rep = factor(rep)) nae_p <- nae_df %>% ggplot(aes(time, rel_N, color = pool)) + geom_hline(yintercept = 0, color = "gray70") + geom_point(data = nae_df %>% filter(time == 0) %>% mutate(time = -5)) + geom_line(size = 0.5) + scale_y_continuous("Proportional change in N", breaks = c(-0.5, 0, 0.5)) + scale_x_continuous("Time (days)", breaks = c(0, 50, 100)) + facet_wrap(~ rep, ncol = 4) + scale_color_manual(NULL, values = color_pal()[c(4:6, 1:3)]) + theme(panel.spacing = unit(1.5, "lines"), strip.text = element_blank()) + NULL nae_p
\clearpage
supp_sims1 <- crossing(.aDV = c(par_estimates$aDV[1], par_estimates$aPH[1]), .aPH = c(par_estimates$aPH[1], par_estimates$aDV[1])) %>% pmap_dfr(function(.aDV, .aPH) { food_web(tmax = 100, s = 10, b = 20, w = 20, other_pars = list(q = 3, aDV = .aDV, aPH = .aPH)) %>% mutate(aDV = .aDV, aPH = .aPH) }) %>% spread(pool, N) %>% mutate(Vg = V_gain(detritivore, detritus, aDV = aDV), Hg = H_gain(plant, herbivore, aPH = aPH)) %>% select(-soil:-midge) %>% gather("variable", "value", Vg:Hg) %>% mutate(pool = factor(ifelse(grepl("^V", variable), "detritivore", "herbivore")), aDV = factor(aDV, levels = range(as.numeric(paste(aDV))), labels = paste(c("low", "high"), "detritivore uptake")), aPH = factor(aPH, levels = range(as.numeric(paste(aPH))), labels = paste(c("low", "high"), "herbivore uptake"))) %>% select(-variable) %>% select(aDV, aPH, pool, everything()) %>% group_by(aDV, aPH, pool) %>% mutate(value = value - value[1]) %>% ungroup()
s1_caption <- "Time series of the bottom-up effects on detritivore and herbivore pools with (a,b) low and (c,d) high herbivore uptake rates ($a_P$), and with (a,c) low and (b,d) high detritivore uptake rates ($a_D$). Parameter values are as in Table 1, with low $a_D = 0.012$, low $a_P = 0.012$, high $a_D = 0.032$, high $a_P = 0.032$, $w = 20$, $b = 20$, and $q = 3$." s1_caption_short <- "Effects of varying herbivore and detritivore uptake rates on bottom-up effects"
trans_p3 <- ggplot(data = NULL) + geom_hline(yintercept = 0, color = "gray70") + geom_line(data = supp_sims1, aes(time, value, color = pool), size = 1, linetype = 1) + scale_y_continuous(expression("Bottom-up effect on pool (" * day^{-1} * ")" )) + scale_x_continuous("Time (days)") + scale_color_manual(values = color_pal()[1:2]) + geom_text(data = tibble(time = c( 35, 28), value = c(0.070, 0.007), aDV = sort(unique(supp_sims1$aDV))[1], aPH = sort(unique(supp_sims1$aPH))[1], lab = c("detritivore", "herbivore")), aes(time, value, label = lab), color = color_pal()[1:2], hjust = 0, vjust = 0, lineheight = 0.75, size = 10 / 2.835, parse = TRUE) + facet_grid(aPH ~ aDV) + theme(legend.position = "none", strip.text = element_text(face = "plain", size = 11, margin = margin(b = 4, l = 4)), panel.spacing.x = unit(2.5, "lines"), strip.background = element_blank(), axis.text = element_text(size = 10, color = "black"), axis.title = element_text(size = 12)) + NULL trans_p3
vain_equil_caption <- "Equilibrium N pool size for all six pools versus the inorganic nutrient input rate ($i_I$). Line type indicates the density-independent loss of soil nutrients ($\\mu_I$). Parameter values are as in Table 1." vain_equil_caption_short <- "Effects of the inorganic nutrient input rate and density-independent loss of soil nutrients on equilibrium N pool sizes"
vary_iI_muI_equil <- function(.) { Z <- equil_pools(iI = .$iI[1], muI = .$muI[1]) mutate(Z[["vals"]], iI = .$iI[1], muI = .$muI[1]) } # Takes ~12 sec with 4 threads vain_equil_df <- crossing(iI = seq(0, 20, length.out = 61), muI = c(0.005, 0.01, 0.02)) %>% split(1:nrow(.)) %>% mclapply(vary_iI_muI_equil, mc.cores = parallel::detectCores()) %>% bind_rows() # This just sets the bounds for the facets: vain_dummy_equil <- tibble(iI = 20, pool = factor(c("soil", "detritus", "plant", upper_levels), levels = c("soil", "detritus", "plant", upper_levels)), eq = c(700, 60, 60, 2, 2, 2)) vain_equil_p <- vain_equil_df %>% mutate_at(vars(muI), factor) %>% ggplot(aes(iI, eq)) + geom_hline(yintercept = 0, color = "gray70") + geom_line(aes(color = pool, linetype = muI), size = 0.75) + geom_blank(data = vain_dummy_equil) + scale_y_continuous(expression("Equilibrium pool size (g N" ~ m^{-2} * ")")) + scale_x_continuous(expression(paste("Inorganic nutrient input rate", ~ "(" * g ~ N ~ m^{-2} ~ d^{-1} * ")"))) + facet_wrap(~ pool, nrow = 2, scales = "free_y") + scale_color_manual(NULL, values = color_pal()[c(4:6, 1:3)], guide = FALSE) + scale_linetype_manual("Density-independent\nloss of soil nutrients:", values = c(2, 1, 3)) + theme(panel.spacing = unit(1, "lines"), legend.position = "top") + NULL vain_equil_p
vain_ts_caption <- "Time series of proportional changes in N content among all six pools with varying values of $i_I$ and $\\mu_I$, the inorganic nutrient input rate and density-independent loss of soil nutrients, respectively. Sub-panel labels indicate the values of each parameter, and middle values are the defaults used in the main text. Lines indicate the N content at time $t$ relative to the starting value [i.e. $(N_t - N_{0}) / N_{0}$]. Parameter values are as in Table 1, with predator exploitation $q = 3$, pulse duration $w = 0$, and pulse rate $b = 0$. Note that when $i_I = 5$ and $\\mu_I = 0.02$, predators go extinct so are not shown in that panel." vain_ts_caption_short <- "Time series of N content in response to a midge pulse, while varying the inorganic nutrient input rate and density-independent loss of soil nutrients"
vary_iI_muI <- function(.iI, .muI) { # .iI <- 5; .muI <- 0.02 # rm(.iI, .muI, pools, sim) if (.iI == 5 && .muI == 0.02) { pools <- equil_pools(iI = .iI, muI = .muI, Xeq = 0) sim <- food_web(tmax = 100, s = 10, b = 50, w = 20, other_pars = list(q = 3, iI = .iI, muI = .muI), ep_obj = pools) %>% filter(pool != "midge") %>% mutate(pool = droplevels(pool)) %>% filter(pool != "predator") } else { pools <- equil_pools(iI = .iI, muI = .muI) sim <- food_web(tmax = 100, s = 10, b = 50, w = 20, other_pars = list(q = 3, iI = .iI, muI = .muI), ep_obj = pools) %>% filter(pool != "midge") %>% mutate(pool = droplevels(pool)) } sim <- sim %>% group_by(pool) %>% mutate(N_rel = (N - N[1]) / (N[1])) %>% ungroup() %>% mutate(iI = .iI, muI = .muI) return(sim) } vain_df <- crossing(.iI = c(5, 10, 20), .muI = c(0.005, 0.01, 0.02)) %>% pmap_dfr(vary_iI_muI) %>% mutate_at(vars(iI, muI), factor) %>% mutate(muI = factor(muI, levels = levels(muI), labels = paste("italic(mu[I]) ==", levels(muI))), iI = factor(iI, levels = levels(iI), labels = paste("italic(i[I]) ==", levels(iI)))) vain_dummy <- vain_df %>% group_by(iI, muI) %>% filter(N_rel == max(N_rel)) %>% ungroup() %>% select(-pool, -N) %>% mutate(time = 10, time2 = 10 + 20, N_rel = -0.1 * N_rel) vain_p <- vain_df %>% ggplot(aes(time, N_rel)) + geom_hline(yintercept = 0, color = "gray70") + geom_segment(data = vain_dummy, aes(xend = time2, yend = N_rel), size = 1.5) + geom_line(aes(color = pool), size = 0.5) + scale_y_continuous("Proportional change in N") + scale_x_continuous("Time (days)") + facet_wrap(iI ~ muI, nrow = 3, scales = "free_y", labeller = label_parsed) + scale_color_manual(NULL, values = color_pal()[c(4:6, 1:3)]) + guides(color = guide_legend(override.aes = list(size = 1))) + theme(panel.spacing = unit(1, "lines"), strip.text = element_text(family = "CMU Serif")) + NULL vain_p
fig2_abs_caption <- "Time series of N content among all six pools. Lines indicate the N content at time $t$. The lower black bar represents the duration of the midge N pulse. Parameter values are as in Table 1, with predator exploitation $q = 3$, pulse duration $w = 20$, and pulse rate $b = 50$." fig2_abs_caption_short <- "Time series of absolute N content in response to a midge pulse"
middle_sim <- food_web(tmax = 100, s = 10, b = 50, w = 20, other_pars = list(q = 3)) # This just sets the bounds for the facets: dummy_data <- tibble(time = 20, pool = factor(rep(c("soil", "detritus", "plant", upper_levels), each = 2), levels = c("soil", "detritus", "plant", upper_levels)), N = c(middle_sim %>% filter(pool == "soil") %>% .[["N"]] %>% max() %>% rep(2), middle_sim %>% filter(pool %in% c("detritus", "plant")) %>% .[["N"]] %>% max() %>% rep(2 * 2), middle_sim %>% filter(pool %in% upper_levels) %>% .[["N"]] %>% max() %>% rep(3 * 2))) %>% group_by(pool) %>% mutate(N = N[1] * c(-0.1, 1)) %>% ungroup() trans_pS1 <- middle_sim %>% filter(pool != "midge") %>% mutate(pool = droplevels(pool)) %>% ggplot(aes(time, N)) + geom_hline(yintercept = 0, color = "gray70") + geom_segment(data = tibble(time = 10, time2 = 10+20, N = 0), aes(xend = time2, yend = N), size = 1.5) + geom_line(aes(color = pool), size = 1) + geom_blank(data = dummy_data) + scale_y_continuous(expression("Nitrogen" ~ (g ~ N ~ m^{-2}))) + scale_x_continuous("Time (days)") + facet_wrap(~ pool, nrow = 2, scales = "free") + scale_color_manual(values = color_pal()[c(4:6, 1:3)], guide = FALSE) + theme(panel.spacing = unit(1.5, "lines")) + geom_text(data = dummy_data %>% filter(N < 0, pool == "soil"), label = "pulse", size = 10 / 2.835, hjust = 0.5, vjust = 0, color = "black") + NULL trans_pS1
\clearpage
nohv_fig3_caption <- "Time series of proportional changes in N content for detritivore and herbivore pools (indicated by the left y-axis). This is shown for when (a) no detritivores are present, (b) both are present, and (c) no herbivores are present. (a,c) When either detritivores or herbivores were absent, we solved for the other equilibrium pool sizes before running the simulations. The blue shaded regions represent the proportional change in N for predators, which is indicated by the right y-axis. Different y-axis scales are used to aid visualization of the transient dynamics of the herbivores and detritivores, which had a weaker response to the pulse than predators for the selected parameter values. The lower black bar represents the duration of the midge N pulse. Parameter values are as in Table 1, with predator exploitation $q = 3$, pulse duration $w = 20$, and pulse rate $b = 20$." nohv_fig3_caption_short <- "Time series of N content in response to a midge pulse, with herbivores or detritivores removed"
do_nohv_fig3 <- function(.H, .V) { stopifnot(length(.H) == 1 && .H %in% 0L:1L) stopifnot(length(.V) == 1 && .V %in% 0L:1L) if (.H == 0 && .V == 0) { pools <- equil_pools(Heq = 0, Veq = 0) } if (.H == 0 && .V == 1) { pools <- equil_pools(Heq = 0) } if (.H == 1 && .V == 0) { pools <- equil_pools(Veq = 0) } if (.H == 1 && .V == 1) { pools <- NULL } q_ <- 3 sim <- food_web(tmax = 100, s = 10, b = 20, w = 20, other_pars = list(q = q_), ep_obj = pools) %>% filter(pool %in% c(upper_levels, "midge")) %>% mutate(pool = droplevels(pool)) %>% mutate(pool = factor(paste(pool), levels = c(upper_levels, "midge")), H = .H, V = .V) return(sim) } nohv_fig3_df <- crossing(.H = 0L:1L, .V = 0L:1L) %>% filter(!(.H == 0 & .V == 0)) %>% pmap_dfr(do_nohv_fig3) %>% group_by(H, V, pool) %>% mutate(N_rel = (N - N[1]) / N[1]) %>% mutate(N_rel = ifelse(is.nan(N_rel), 0, N_rel)) %>% ungroup() nohv_fig3_plot <- function(.H, .V, .mult = 25.2, .ylims = c(-0.3011, 1.1)) { # .H = 1; .V = 0 # .mult = 25.2; .ylims = c(-0.3011, 1.1) # rm(.H, .V, .mult, .ylims) .title <- case_when(.H == 0 & .V == 1 ~ "No herbivores", .H == 1 & .V == 0 ~ "No detritivores", .H == 1 & .V == 1 ~ "Both present", TRUE ~ "ERROR") hv_combo <- nohv_fig3_df %>% filter(H == .H, V == .V) %>% select(-H, -V) hv_plot <- hv_combo %>% filter(!pool %in% c("midge", "predator")) %>% mutate(pool = droplevels(pool)) %>% ggplot(aes(time, N_rel)) + geom_hline(yintercept = 0, color = "gray70") + geom_segment(data = tibble(time = 10, time2 = 10+20, N_rel = -0.1), aes(xend = time2, yend = N_rel), size = 1.5) + geom_ribbon(data = hv_combo %>% filter(pool == "predator") %>% mutate(N_rel = N_rel / .mult), aes(ymin = 0, ymax = N_rel), fill = color_pal(0.5)[3]) + geom_line(aes(color = pool), size = 0.75) + scale_x_continuous("Time (days)") + scale_color_manual(NULL, values = color_pal()) + ggtitle(.title) + theme(legend.position = "none", strip.text.y = element_text(face = "plain", size = 10, angle = 270, margin = margin(l = 4)), strip.text.x = element_text(face = "plain", size = 10, margin = margin(b = 4)), panel.spacing.y = unit(0, "lines"), panel.spacing.x = unit(1.5, "lines"), plot.margin = margin(0,0,t=12,b=12), plot.title = element_text(size = 12, hjust = 0.5, face = "plain", margin = margin(0,0,0,b = 10))) + scale_y_continuous("Proportional change in N", limits = .ylims, breaks = c(0, 0.5, 1), sec.axis = sec_axis(~ . * .mult, name = paste("Proportional change in N", "(predator)"), breaks = c(0, 10, 20))) + NULL if (.H == 1 && .V == 1) { hv_plot <- hv_plot + geom_text(data = tibble(pool = factor(upper_levels[upper_levels != "predator"]), time = c(92, 100), N_rel = c(0.3, -0.1), q = factor(paste(rep("high", 2), "exploitation"), levels = paste(c("low", "high"), "exploitation"))), aes(label = pool, color = pool), hjust = 1, size = 10 / 2.835) + geom_text(data = tibble(time = 20, N_rel = -0.15, q = factor("high exploitation", levels = paste(c("low", "high"), "exploitation"))), label = "pulse", size = 10 / 2.835, hjust = 0.5, vjust = 1, color = "black") + theme(axis.title.y.left = element_text(margin = margin(0,0,0,r=12), size = 11), axis.title.y.right = element_text(margin = margin(0,0,0,l=12), size = 11)) } else { hv_plot <- hv_plot + theme(axis.title.y = element_blank()) } if (.H == 1 && .V == 0) { hv_plot <- hv_plot + geom_text(data = tibble(time = 65, N_rel = 0.8, q = factor("high exploitation", levels = paste(c("low", "high"), "exploitation"))), label = "predator", color = color_pal()[3], hjust = 1, vjust = 1, size = 10 / 2.835) } if (!(.H == 0 & .V == 1)) { hv_plot <- hv_plot + theme(axis.title.x = element_blank(), axis.text.x = element_blank()) } else { hv_plot <- hv_plot + theme(axis.title.x = element_text(margin = margin(0,0,0,t=4))) } return(hv_plot) } nohv_fig3_plot_list <- crossing(.H = 0L:1L, .V = 0L:1L) %>% filter(!(.H == 0 & .V == 0)) %>% .[c(2, 3, 1),] %>% pmap(nohv_fig3_plot) %>% map(~ .x + theme(axis.line.y.left = element_line(size = 1))) ggarrange(plots = nohv_fig3_plot_list, ncol = 1, labels = sprintf("(%s)", letters[1:3]), label.args = list(gp = gpar(fontsize = 14, fontface = "plain"), vjust = 2, hjust = -1.5))
nohv_fig4_caption <- "Time series of the top-down (\"TD\") and bottom-up (\"BU\") effects on detritivore and herbivore pools. This is shown for when (a) no detritivores are present, (b) both are present, and (c) no herbivores are present. (a,c) When either detritivores or herbivores were absent, we solved for the other equilibrium pool sizes before running the simulations. The gray line indicates the top-down effects on both detritivore and herbivore pools, as they are identical in the model. Parameter values are as in Table 1, with predator exploitation $q = 3$, pulse duration $w = 20$, and pulse rate $b = 20$." nohv_fig4_caption_short <- "Time series of top-down and bottom-up effects, with herbivores or detritivores removed"
do_nohv_fig4 <- function(.H, .V) { stopifnot(length(.H) == 1 && .H %in% 0L:1L) stopifnot(length(.V) == 1 && .V %in% 0L:1L) if (.H == 0 && .V == 0) { pools <- equil_pools(Heq = 0, Veq = 0) } if (.H == 0 && .V == 1) { pools <- equil_pools(Heq = 0) } if (.H == 1 && .V == 0) { pools <- equil_pools(Veq = 0) } if (.H == 1 && .V == 1) { pools <- NULL } q_ <- 3 sim <- food_web(tmax = 100, s = 10, b = 20, w = 20, other_pars = list(q = q_), ep_obj = pools) %>% spread(pool, N) %>% mutate(Vg = V_gain(detritivore, detritus), Vl = V_loss(detritivore, predator, herbivore, midge, q_), Hg = H_gain(plant, herbivore), Hl = H_loss(herbivore, predator, detritivore, midge, q_)) %>% select(-soil:-midge) %>% gather("variable", "value", Vg:Hl) %>% mutate(pool = factor(ifelse(grepl("^V", variable), "detritivore", "herbivore")), type = factor(ifelse(grepl("l$", variable), "top-down", "bottom-up"))) %>% select(-variable) %>% select(pool, type, everything()) %>% group_by(pool, type) %>% mutate(value = value - value[1]) %>% ungroup() %>% mutate(H = .H, V = .V) return(sim) } nohv_fig4_df <- crossing(.H = 0L:1L, .V = 0L:1L) %>% filter(!(.H == 0 & .V == 0)) %>% pmap_dfr(do_nohv_fig4) nohv_fig4_plot <- function(.H, .V, .ylims = c(-0.037, 0.144)) { .title <- case_when(.H == 0 & .V == 1 ~ "No herbivores", .H == 1 & .V == 0 ~ "No detritivores", .H == 1 & .V == 1 ~ "Both present", TRUE ~ "ERROR") .td_pool <- ifelse(.V == 0, "herbivore", "detritivore") hv_combo <- nohv_fig4_df %>% filter(H == .H, V == .V) %>% select(-H, -V) %>% mutate(value = ifelse(is.nan(value), 0, value)) hv_plot <- ggplot(data = NULL) + geom_hline(yintercept = 0, color = "gray70") + geom_line(data = hv_combo %>% filter(type == "bottom-up"), aes(time, value, color = pool, group = interaction(pool, type)), size = 1, linetype = 1) + geom_line(data = hv_combo %>% filter(type == "top-down", pool == .td_pool), aes(time, value), size = 1, color = "gray60") + xlab("Time (days)") + ggtitle(.title) + scale_color_manual(values = color_pal()[1:2]) + scale_y_continuous(expression("Effect on pool (" * day^{-1} * ")" ), breaks = c(0, 0.06, 0.12), limits = .ylims) + theme(legend.position = "none", strip.text = element_text(face = "plain", size = 11, margin = margin(b = 4)), panel.spacing.x = unit(1.5, "lines"), plot.margin = margin(0,t=12,b=12,r=4), strip.background = element_blank(), axis.text = element_text(size = 10, color = "black"), axis.title = element_text(size = 12), plot.title = element_text(size = 12, hjust = 0.5, face = "plain", margin = margin(0,0,0,b=6))) if (.H == 1 && .V == 1) { hv_plot <- hv_plot + geom_text(data = tibble(time = c( 50), value = c(0.05), lab = sprintf("italic(%s)", c("'BU'['detritivore']"))), aes(time, value, label = lab), color = c(color_pal()[1]), hjust = 0, vjust = 0, lineheight = 0.75, size = 10 / 2.835, parse = TRUE) + theme(axis.title.y = element_text(margin = margin(0,0,r=4,l=10))) } else { hv_plot <- hv_plot + theme(axis.title.y = element_blank()) } if (.H == 1 && .V == 0) { hv_plot <- hv_plot + geom_text(data = tibble(time = c(5, 50), value = c(0.034, 0.05), lab = sprintf("italic(%s)", c("'BU'['herbivore']", "'TD'['both']"))), aes(time, value, label = lab), color = c(color_pal()[2], "gray60"), hjust = 0, vjust = 0, lineheight = 0.75, size = 10 / 2.835, parse = TRUE) } if (!(.H == 0 & .V == 1)) { hv_plot <- hv_plot + theme(axis.title.x = element_blank(), axis.text.x = element_blank()) } else { hv_plot <- hv_plot + theme(axis.title.x = element_text(margin = margin(0,0,0,t=4))) } return(hv_plot) } nohv_fig4_plot_list <- crossing(.H = 0L:1L, .V = 0L:1L) %>% filter(!(.H == 0 & .V == 0)) %>% .[c(2, 3, 1),] %>% pmap(nohv_fig4_plot) %>% map(~ .x + theme(axis.line.y.left = element_line(size = 1))) ggarrange(plots = nohv_fig4_plot_list, ncol = 1, labels = sprintf("(%s)", letters[1:3]), label.args = list(gp = gpar(fontsize = 14, fontface = "plain"), vjust = 2, hjust = -2))
mtop_fig5_caption <- "Cumulative midge effect on total bottom-up control on detritivores (${}^V BU_{total}$) total bottom-up control on herbivores (${}^H BU_{total}$), net top-down control ($TD_{net}$), top-down intensification, and top-down alleviation as a function of total midge inputs. This is shown for when (a) midges only go to the detritus pool, (b) midges go to both detritus and predator pools, and (c) midges only go to the predator pool. The panels show different combinations of high and low predator exploitation and handling times on midges. Parameter values are as in Table 1, with low exploitation $q = 0.8$, high exploitation $q = 8$, low midge handling time $h_M = 1.3$, high midge handling time $h_M = 5.3$, pulse duration $w = 20$, and pulse rate $b$ ranging from 0.1 to 23 to produce different total midge inputs." mtop_fig5_caption_short <- "Total effect of midges on bottom-up and top-down control, with varying midge exploitation and handling time, and with midges going only to detritus or to predators"
do_mtop_fig5 <- function(.midges_not_to) { # .midges_not_to = "none" # rm(.midges_not_to) stopifnot(length(.midges_not_to) == 1 && .midges_not_to %in% c("X", "none", "D")) .one_combo <- function(.q, .hM, .b) { # .q = 0.8; .hM = par_estimates$hM[1] * 0.5; .b = 0.1 # rm(.q, .hM, .b) fw <- food_web(tmax = 500, s = 10, b = .b, w = 20, other_pars = list(hM = .hM, q = .q), midges_not_to = .midges_not_to) fw <- fw %>% spread(pool, N) %>% mutate(Vg = V_gain(detritivore, detritus), Vl = V_loss(detritivore, predator, herbivore, midge, q = .q, hM = .hM, .no_M_to_X = (.midges_not_to == "X")), Hg = H_gain(plant, herbivore), Hl = H_loss(herbivore, predator, detritivore, midge, q = .q, hM = .hM, .no_M_to_X = (.midges_not_to == "X"))) %>% summarize(cum_pos_loss_V = sum(Vl[Vl > Vl[1]] - Vl[1]), cum_pos_loss_H = sum(Hl[Hl > Hl[1]] - Hl[1]), cum_neg_loss_V = sum(Vl[Vl < Vl[1]] - Vl[1]), cum_neg_loss_H = sum(Hl[Hl < Hl[1]] - Hl[1]), cum_gain_V = sum(Vg[Vg > Vg[1]] - Vg[1]), cum_gain_H = sum(Hg[Hg > Hg[1]] - Hg[1])) %>% mutate(area = 20 * .b, q = .q, hM = .hM) %>% select(area, q, hM, everything()) %>% mutate(cum_loss_V = cum_pos_loss_V + cum_neg_loss_V, cum_loss_H = cum_pos_loss_H + cum_neg_loss_H, # Changing to magnitudes cum_neg_loss_V = abs(cum_neg_loss_V), cum_neg_loss_H = abs(cum_neg_loss_H)) %>% gather("variable", "value", cum_pos_loss_V:cum_loss_H) %>% mutate(pool = factor(ifelse(grepl("V$", variable), "detritivore", "herbivore")), direction = factor(ifelse(grepl("loss", variable), "top-down", "bottom-up")), type = case_when(grepl("^cum_pos_loss", variable) ~ "TD intensification", grepl("^cum_neg_loss", variable) ~ "TD alleviation", grepl("^cum_loss_", variable) ~ "TD net", grepl("^cum_gain_", variable) ~ "BU total", TRUE ~ NA_character_) %>% factor(levels = c("BU total", "TD net", "TD intensification", "TD alleviation"))) %>% select(-variable) if (sum(is.na(fw$type)) > 0) stop("Error in processing `type` column.") return(fw) } out <- crossing(.q = c(0.8, 8), .hM = par_estimates$hM[1] * c(0.5, 2), .b = seq(0.1, 23, length.out = 21)) %>% pmap_dfr(.one_combo) %>% mutate(q = factor(q, levels = sort(unique(q)), labels = sprintf("%s midge\nexploitation", c("low", "high"))), hM = factor(hM, levels = sort(unique(hM)), labels = paste(c("low", "high"), "midge\nhandling time"))) %>% mutate(not_to = .midges_not_to) return(out) } # Takes ~24 sec using 3 threads mtop_fig5_df <- mclapply(c("X", "none", "D"), do_mtop_fig5, mc.cores = min(3, parallel::detectCores())) %>% bind_rows() mtop_fig5_plot <- function(.midges_not_to, .ylims = c(-0.688, 4.305)) { fig5_labs <- list(bquote(italic({}^'V' * 'BU'['total'])), bquote(italic({}^'H' * 'BU'['total'])), bquote(italic('TD'['net'])), bquote(italic('TD'['intensification'])), bquote(italic('TD'['alleviation']))) fig5_lvls <- c("detritivore BU total", "herbivore BU total", "TD net", "TD intensification", "TD alleviation") .colors <- c(color_pal()[1:2], "gray60", "gray30", "gray60") # .midges_not_to = "D"; .ylims = c(-0.688, 4.305) # rm(.midges_not_to, .ylims, .pool) .title <- case_when(.midges_not_to == "none" ~ "Midges to detritus and predators", .midges_not_to == "D" ~ "Midges to predators only", .midges_not_to == "X" ~ "Midges to detritus only", TRUE ~ "ERROR") mt_combo <- mtop_fig5_df %>% filter(not_to == .midges_not_to) %>% select(-not_to) %>% mutate(id = ifelse(grepl("^TD", type), paste(type), sprintf("%s %s", paste(pool), paste(type))), id = factor(id, levels = fig5_lvls)) %>% group_by(id, area, q, hM) %>% summarize(value = value[1]) %>% ungroup() mt_plot <- mt_combo %>% ggplot(aes(area, value)) + geom_hline(yintercept = 0, color = "gray80") + geom_line(aes(color = id, linetype = id), size = 0.75) + ggtitle(.title) + facet_grid(hM ~ q) + scale_color_manual(NULL, values = .colors, guide = FALSE) + scale_linetype_manual(NULL, values = c(1, 1, 1, 3:2), breaks = fig5_lvls, labels = fig5_labs) + guides(linetype = guide_legend(keywidth = 2, ncol = 1, keyheight = 0.6, override.aes = list(color = .colors, size = 1), byrow = TRUE)) + xlab(expression("Total midge input (" * g ~ N ~ m^{-2} * ")")) + scale_y_continuous("Total midge effect", limits = .ylims, breaks = c(0, 2, 4)) + theme(legend.background = element_blank(), legend.position = "right", legend.direction = "vertical", legend.justification = "left", legend.box = "vertical", legend.title = element_text(size = 11, margin = margin(0,0,0,b=-4)), legend.text = element_text(size = 10, family = "CMU Serif"), strip.background = element_blank(), strip.text.y = element_text(vjust = 0.5, hjust = 0, size = 9, angle = 0), strip.text.x = element_text(size = 9), panel.spacing.x = unit(1.5, "lines"), plot.margin = margin(0,t=12,b=12,r=4), plot.title = element_text(size = 12, hjust = 0.5, face = "plain", margin = margin(0,0,0,b=6))) + NULL if (.midges_not_to == "none") { mt_plot <- mt_plot + theme(axis.title.y = element_text(margin = margin(0,0,0,r=12), size = 12)) } else { mt_plot <- mt_plot + theme(axis.title.y = element_blank(), legend.position = "none") } if (.midges_not_to != "D") { mt_plot <- mt_plot + theme(axis.title.x = element_blank(), axis.text.x = element_blank()) } else { mt_plot <- mt_plot + theme(axis.title.x = element_text(margin = margin(0,0,0,t=4), size = 12)) } return(mt_plot) } mtop_fig5_plot_list <- tibble(.midges_not_to = c("X", "none", "D")) %>% pmap(mtop_fig5_plot) ggarrange(plots = mtop_fig5_plot_list, ncol = 1, labels = sprintf("(%s)", letters[1:3]), label.args = list(gp = gpar(fontsize = 14, fontface = "plain"), vjust = 2, hjust = -1))
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.