This script presents R code for the analyses and figures of our study Divergent roles of herbivory in eutrophying forests. Please see the manuscript for details on the rationale and explanation of these analyses.
For all models, the following abbreviations apply:
# load packages
library(tidyverse)
library(tidybayes)
library(brms)
library(patchwork)
library(hrbrthemes)
library(see)
# source data 1
d1 <- read_csv("data/source_data1.csv")
d1
## # A tibble: 52 × 8
## id d_herbivory baseline_herbiv…¹ time_…² site_…³ herb_…⁴ shrub…⁵ tree_…⁶
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 EU_006 0 6 46 3.12 -66.7 -11.9 50.4
## 2 EU_009a 6 1 41 1.56 -14.1 -17.9 2.28
## 3 EU_009b 6 1 26 1.56 -47.0 -3.90 -7.50
## 4 EU_011 0 6 22 5.88 -78.8 -12.7 -32.8
## 5 EU_012 6 11 49 2.54 -56.3 0.521 19.3
## 6 EU_013 18 7 53 3.32 -36.4 -19.7 -2.26
## 7 EU_014 6 6 56 3.68 -94.0 1.44 -16.9
## 8 EU_015 11 5 25 2.53 24.9 NA -15.6
## 9 EU_016 0 2 21 1.08 93.6 18.0 2.62
## 10 EU_025 2 7 51 3.60 46.5 26.4 2.06
## # … with 42 more rows, and abbreviated variable names ¹baseline_herbivory,
## # ²time_span, ³site_area_log, ⁴herb_cover, ⁵shrub_cover, ⁶tree_cover
## # ℹ Use `print(n = ...)` to see more rows
# source data 2
d2 <- read_csv("data/source_data2.csv")
d2
## # A tibble: 52 × 7
## id d_herbivory baseline_herbivory time_span site_area_…¹ d_spp…² spp_e…³
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 EU_006 0 6 46 3.12 -4 0.505
## 2 EU_009a 6 1 41 1.56 -21 0.536
## 3 EU_009b 6 1 26 1.56 5 0.593
## 4 EU_011 0 6 22 5.88 26 0.263
## 5 EU_012 6 11 49 2.54 -93 0.546
## 6 EU_013 18 7 53 3.32 -32 0.524
## 7 EU_014 6 6 56 3.68 -33 0.488
## 8 EU_015 11 5 25 2.53 -18 0.356
## 9 EU_016 0 2 21 1.08 -2 0.3
## 10 EU_025 2 7 51 3.60 -15 0.5
## # … with 42 more rows, and abbreviated variable names ¹site_area_log,
## # ²d_spp_richness, ³spp_exchange_ratio
## # ℹ Use `print(n = ...)` to see more rows
# source data 3
d3 <- read_csv("data/source_data3.csv")
d3
## # A tibble: 52 × 10
## id d_herbivory baseli…¹ time_…² site_…³ d_smal…⁴ d_cwmn d_ali…⁵ d_red…⁶
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 EU_006 0 6 46 3.12 -0.114 0.144 5.29 0
## 2 EU_009a 6 1 41 1.56 1.53 0.184 2.20 5.94
## 3 EU_009b 6 1 26 1.56 0.278 -0.0742 3.06 -3.61
## 4 EU_011 0 6 22 5.88 0.0942 -0.0324 0.541 -0.306
## 5 EU_012 6 11 49 2.54 -9.26 0.587 -0.0987 -3.45
## 6 EU_013 18 7 53 3.32 -19.4 0.888 11.1 -5.85
## 7 EU_014 6 6 56 3.68 -8.16 -0.216 0.590 0.590
## 8 EU_015 11 5 25 2.53 0.249 -0.112 0 1.29
## 9 EU_016 0 2 21 1.08 -4.44 0.0613 0 0.0264
## 10 EU_025 2 7 51 3.60 2.49 0.266 0.898 -3.82
## # … with 42 more rows, 1 more variable: cum_ndepo <dbl>, and abbreviated
## # variable names ¹baseline_herbivory, ²time_span, ³site_area_log,
## # ⁴d_small_ranged_percent, ⁵d_alien_percent, ⁶d_redlist_percent
## # ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable names
mod_shrub <- brms::brm(data = d1,
scale(shrub_cover) ~
scale(d_herbivory) +
scale(baseline_herbivory) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 2000, chains = 4, cores = 4,
file = "models/mod_shrub")
mod_shrub
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(shrub_cover) ~ scale(d_herbivory) + scale(baseline_herbivory) + scale(time_span) + scale(site_area_log)
## Data: d1 (Number of observations: 50)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -0.02 0.13 -0.28 0.25 1.00 5130
## scaled_herbivory -0.42 0.17 -0.74 -0.10 1.00 4811
## scalebaseline_herbivory 0.13 0.15 -0.16 0.43 1.00 4494
## scaletime_span 0.14 0.15 -0.17 0.44 1.00 4469
## scalesite_area_log -0.01 0.16 -0.31 0.30 1.00 5084
## Tail_ESS
## Intercept 2170
## scaled_herbivory 3099
## scalebaseline_herbivory 3192
## scaletime_span 2765
## scalesite_area_log 3124
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.96 0.10 0.79 1.18 1.00 4365 2823
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for negative slope
mod_shrub %>% spread_draws(b_scaled_herbivory) %>% filter(b_scaled_herbivory < 0) %>% nrow /
mod_shrub %>% spread_draws(b_scaled_herbivory) %>% nrow
## [1] 0.9955
ce_mod_shrub <- conditional_effects(mod_shrub, "d_herbivory")
# regression
plot(
ce_mod_shrub,
plot = FALSE,
points = TRUE,
point_args = list(
fill = "#ff00bf",
alpha = 0.6,
pch = 21,
size = 3
),
line_args = list(
col = "black",
fill = "#ff00bf",
alpha = 0.3,
lwd = .6,
lty = 1
)
)[[1]] +
scale_fill_manual(values = c("#ff00bf")) +
scale_colour_manual(values = c("#ff00bf")) +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("\u0394 shrub cover") +
xlab("\u0394 herbivory pressure") -> fig2a1
# boxplot
ggplot(data = d1,
aes(y = shrub_cover, x = 1)) +
geom_violinhalf(fill = "#ff00bf", alpha = 0.4) +
geom_point(
aes(y = shrub_cover, x = 0.92),
position = position_jitter(width = 0.04),
col = "#ff00bf",
alpha = 0.3
) +
geom_boxplot(
aes(x = 0.92) ,
outlier.shape = NA,
width = 0.1,
alpha = 0.4,
fill = "#b631cb"
) +
stat_summary(
fun = mean,
geom = "point",
shape = 24,
size = 6,
color = "black",
fill = "#ff00bf",
alpha = 0.8
) +
geom_hline(yintercept = 0) +
labs(x = "Density", y = "") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.text.x = element_text(colour = 'white'),
axis.text.y = element_text(colour = 'white'),
aspect.ratio = 2
) +
scale_y_continuous(position = "right") -> fig2a2
fig2a1 + fig2a2
mod_herb <- brms::brm(data = d1,
scale(herb_cover) ~
scale(d_herbivory) +
scale(baseline_herbivory) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 2000, chains = 4, cores = 4,
file = "models/mod_herb")
mod_herb
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(herb_cover) ~ scale(d_herbivory) + scale(baseline_herbivory) + scale(time_span) + scale(site_area_log)
## Data: d1 (Number of observations: 51)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 0.01 0.14 -0.28 0.29 1.00 5925
## scaled_herbivory -0.02 0.17 -0.36 0.32 1.00 4494
## scalebaseline_herbivory -0.09 0.16 -0.40 0.22 1.00 5682
## scaletime_span -0.26 0.16 -0.57 0.06 1.00 5240
## scalesite_area_log 0.11 0.17 -0.21 0.44 1.00 5153
## Tail_ESS
## Intercept 3413
## scaled_herbivory 3098
## scalebaseline_herbivory 2645
## scaletime_span 2813
## scalesite_area_log 3062
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.03 0.11 0.84 1.27 1.00 4253 3037
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for negative slope
mod_herb %>% spread_draws(b_scaled_herbivory) %>% filter(b_scaled_herbivory < 0) %>% nrow /
mod_herb %>% spread_draws(b_scaled_herbivory) %>% nrow
## [1] 0.55275
ce_mod_herb <- conditional_effects(mod_herb, "d_herbivory")
# regression
plot(
ce_mod_herb,
plot = FALSE,
points = TRUE,
point_args = list(
fill = "#bfff00",
alpha = 0.6,
pch = 21,
size = 3
),
line_args = list(
col = "black",
fill = "#bfff00",
alpha = 0.3,
lwd = .6,
lty = 2
)
)[[1]] +
scale_fill_manual(values = c("#bfff00")) +
scale_colour_manual(values = c("#bfff00")) +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("\u0394 herb cover") +
xlab("\u0394 herbivory pressure") -> fig2b1
# boxplot
ggplot(data = d1, aes(y = herb_cover, x = 1)) +
geom_violinhalf(fill = "#bfff00", alpha = 0.4) +
geom_point(
aes(y = herb_cover, x = 0.92),
position = position_jitter(width = 0.04),
col = "#bfff00",
alpha = 0.3
) +
geom_boxplot(
aes(x = 0.92) ,
outlier.shape = NA,
width = 0.1,
alpha = 0.4,
fill = "#bfff00"
) +
geom_hline(yintercept = 0) +
stat_summary(
fun = mean,
geom = "point",
shape = 24,
size = 6,
color = "black",
fill = "#bfff00",
alpha = 0.8
) +
labs(x = "Density", y = "") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.text.x = element_text(colour = 'white'),
axis.text.y = element_text(colour = 'white'),
aspect.ratio = 2
) +
scale_y_continuous(position = "right") -> fig2b2
fig2b1 + fig2b2
mod_tree <- brms::brm(data = d1,
scale(tree_cover) ~
scale(d_herbivory) +
scale(baseline_herbivory) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 2000, chains = 4, cores = 4,
file = "models/mod_tree")
mod_tree
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(tree_cover) ~ scale(d_herbivory) + scale(baseline_herbivory) + scale(time_span) + scale(site_area_log)
## Data: d1 (Number of observations: 51)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -0.01 0.15 -0.30 0.28 1.00 6000
## scaled_herbivory 0.12 0.17 -0.22 0.48 1.00 4911
## scalebaseline_herbivory 0.12 0.16 -0.21 0.43 1.00 5528
## scaletime_span 0.08 0.17 -0.24 0.42 1.00 5570
## scalesite_area_log 0.05 0.17 -0.27 0.38 1.00 4804
## Tail_ESS
## Intercept 2545
## scaled_herbivory 3114
## scalebaseline_herbivory 3040
## scaletime_span 2931
## scalesite_area_log 3143
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.04 0.11 0.86 1.30 1.00 4351 2643
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for negative slope
mod_tree %>% spread_draws(b_scaled_herbivory) %>% filter(b_scaled_herbivory < 0) %>% nrow /
mod_tree %>% spread_draws(b_scaled_herbivory) %>% nrow
## [1] 0.234
ce_mod_tree <- conditional_effects(mod_tree, "d_herbivory")
# regression
plot(
ce_mod_tree,
plot = FALSE,
points = TRUE,
point_args = list(
fill = "#00bfff",
alpha = 0.6,
pch = 21,
size = 3
),
line_args = list(
col = "black",
fill = "#00bfff",
alpha = 0.3,
lwd = .6,
lty = 2
)
)[[1]] +
scale_fill_manual(values = c("#00bfff")) +
scale_colour_manual(values = c("#00bfff")) +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("\u0394 tree cover") +
xlab("\u0394 herbivory pressure") -> fig2c1
# boxplot
ggplot(data = d1,
aes(y = tree_cover, x = 1)) +
geom_violinhalf(fill = "#00bfff", alpha = 0.4) +
geom_point(
aes(y = tree_cover, x = 0.92),
position = position_jitter(width = 0.04),
col = "#00bfff",
alpha = 0.3
) +
geom_boxplot(
aes(x = 0.92) ,
outlier.shape = NA,
width = 0.1,
alpha = 0.4,
fill = "#00bfff"
) +
stat_summary(
fun = mean,
geom = "point",
shape = 24,
size = 6,
color = "black",
fill = "#00bfff",
alpha = 0.8
) +
geom_hline(yintercept = 0) +
labs(x = "Density", y = "") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.text.x = element_text(colour = 'white'),
axis.text.y = element_text(colour = 'white'),
aspect.ratio = 2
) +
scale_y_continuous(position = "right") -> fig2c2
fig2c1 + fig2c2
fig2a1 + fig2a2 + fig2b1 + fig2b2 + fig2c1 + fig2c2 +
plot_layout(ncol = 4,
nrow = 2,
widths = c(1, 0.5, 1, 0.5)) +
plot_annotation(tag_levels = list(c("a", "", "b", "", "c", ""))) &
theme(plot.tag = element_text(face = 'bold'))
mod_sr <- brms::brm(data = d2,
scale(d_spp_richness) ~
scale(d_herbivory) +
scale(baseline_herbivory) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 2000, chains = 4, cores = 4,
file = "models/mod_sr")
mod_sr
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(d_spp_richness) ~ scale(d_herbivory) + scale(baseline_herbivory) + scale(time_span) + scale(site_area_log)
## Data: d2 (Number of observations: 52)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 0.00 0.14 -0.27 0.27 1.00 6220
## scaled_herbivory -0.23 0.17 -0.56 0.11 1.00 5125
## scalebaseline_herbivory 0.02 0.15 -0.27 0.31 1.00 5434
## scaletime_span -0.27 0.15 -0.58 0.03 1.00 5038
## scalesite_area_log 0.27 0.16 -0.04 0.60 1.00 4954
## Tail_ESS
## Intercept 2808
## scaled_herbivory 3116
## scalebaseline_herbivory 2565
## scaletime_span 2768
## scalesite_area_log 2768
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.99 0.11 0.81 1.22 1.00 4659 3224
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for negative slope
mod_sr %>% spread_draws(b_scaled_herbivory) %>% filter(b_scaled_herbivory < 0) %>% nrow /
mod_sr %>% spread_draws(b_scaled_herbivory) %>% nrow
## [1] 0.911
ce_mod_sr <- conditional_effects(mod_sr, "d_herbivory")
# regression
plot(
ce_mod_sr,
plot = FALSE,
points = TRUE,
point_args = list(
fill = "#ffc000",
alpha = 0.6,
pch = 21,
size = 3
),
line_args = list(
col = "black",
fill = "#ffc000",
alpha = 0.3,
lwd = .6,
lty = 2
)
)[[1]] +
scale_fill_manual(values = c("#ffc000")) +
scale_colour_manual(values = c("#ffc000")) +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("\u0394 species richness") +
xlab("\u0394 herbivory pressure") -> fig3a1
# boxplot
ggplot(data = d2,
aes(y = d_spp_richness, x = 1)) +
geom_violinhalf(fill = "#ffc000", alpha = 0.4) +
geom_point(
aes(y = d_spp_richness, x = 0.92),
position = position_jitter(width = 0.04),
col = "#ffc000",
alpha = 0.3
) +
geom_boxplot(
aes(x = 0.92) ,
outlier.shape = NA,
width = 0.1,
alpha = 0.4,
fill = "#ffc000"
) +
geom_hline(yintercept = 0) +
stat_summary(
fun = mean,
geom = "point",
shape = 24,
size = 6,
color = "black",
fill = "#ffc000",
alpha = 0.8
) +
labs(x = "Density", y = "") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.text.x = element_text(colour = 'white'),
axis.text.y = element_text(colour = 'white'),
aspect.ratio = 2
) +
scale_y_continuous(position = "right") -> fig3a2
fig3a1 + fig3a2
mod_ex <- brms::brm(data = d2,
scale(spp_exchange_ratio) ~
scale(d_herbivory) +
scale(baseline_herbivory) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 2000, chains = 4, cores = 4,
file = "models/mod_ex")
mod_ex
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(spp_exchange_ratio) ~ scale(d_herbivory) + scale(baseline_herbivory) + scale(time_span) + scale(site_area_log)
## Data: d2 (Number of observations: 52)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 0.00 0.12 -0.24 0.24 1.00 5460
## scaled_herbivory 0.52 0.15 0.24 0.81 1.00 4211
## scalebaseline_herbivory 0.04 0.12 -0.20 0.28 1.00 5818
## scaletime_span 0.32 0.13 0.05 0.58 1.00 4942
## scalesite_area_log -0.25 0.14 -0.52 0.03 1.00 5031
## Tail_ESS
## Intercept 2639
## scaled_herbivory 2973
## scalebaseline_herbivory 2994
## scaletime_span 3158
## scalesite_area_log 2847
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.84 0.09 0.69 1.03 1.00 4811 3061
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for positive slope
mod_ex %>% spread_draws(b_scaled_herbivory) %>% filter(b_scaled_herbivory > 0) %>% nrow /
mod_ex %>% spread_draws(b_scaled_herbivory) %>% nrow
## [1] 0.99925
ce_mod_ex <- conditional_effects(mod_ex, "d_herbivory")
# regression
plot(
ce_mod_ex,
plot = FALSE,
points = TRUE,
point_args = list(
fill = "#ff0040",
alpha = 0.6,
pch = 21,
size = 3
),
line_args = list(
col = "black",
fill = "#ff0040",
alpha = 0.3,
lwd = .6,
lty = 1
)
)[[1]] +
scale_fill_manual(values = c("#ff0040")) +
scale_colour_manual(values = c("#ff0040")) +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("Species exchange ratio") +
xlab("\u0394 herbivory pressure") -> fig3b1
# boxplot
ggplot(data = d2,
aes(y = spp_exchange_ratio, x = 1)) +
stat_summary(
fun = mean,
geom = "point",
shape = 24,
size = 6,
color = "black",
fill = "#ff0040",
alpha = 0.8
) +
geom_violinhalf(fill = "#ff0040", alpha = 0.4) +
geom_point(
aes(y = spp_exchange_ratio, x = 0.92),
position = position_jitter(width = 0.04),
col = "#ff0040",
alpha = 0.3
) +
geom_boxplot(
aes(x = 0.92) ,
outlier.shape = NA,
width = 0.1,
alpha = 0.4,
fill = "#ff0040"
) +
stat_summary(
fun = mean,
geom = "point",
shape = 24,
size = 6,
color = "black",
fill = "#ff0040",
alpha = 0.8
) +
labs(x = "Density", y = "") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.text.x = element_text(colour = 'white'),
axis.text.y = element_text(colour = 'white'),
aspect.ratio = 2
) +
scale_y_continuous(position = "right") -> fig3b2
fig3b1 + fig3b2
fig3a1 + fig3a2 + fig3b1 + fig3b2 +
plot_layout(ncol = 4,
widths = c(1, 0.5, 1, 0.5)) +
plot_annotation(tag_levels = list(c("a", "", "b", ""))) &
theme(plot.tag = element_text(face = 'bold'))
mod_cwmn <- brms::brm(data = d3,
scale(d_cwmn) ~
scale(d_herbivory) +
scale(baseline_herbivory) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 2000, chains = 4, cores = 4,
file = "models/mod_cwmn")
mod_cwmn
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(d_cwmn) ~ scale(d_herbivory) + scale(baseline_herbivory) + scale(time_span) + scale(site_area_log)
## Data: d3 (Number of observations: 52)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -0.00 0.12 -0.24 0.24 1.00 5088
## scaled_herbivory 0.43 0.15 0.13 0.73 1.00 3977
## scalebaseline_herbivory 0.01 0.13 -0.24 0.25 1.00 4786
## scaletime_span 0.16 0.14 -0.10 0.43 1.00 4785
## scalesite_area_log 0.12 0.14 -0.16 0.40 1.00 4738
## Tail_ESS
## Intercept 2691
## scaled_herbivory 2836
## scalebaseline_herbivory 2724
## scaletime_span 2626
## scalesite_area_log 3112
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.87 0.10 0.70 1.09 1.00 4498 2735
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for positive slope
mod_cwmn %>% spread_draws(b_scaled_herbivory) %>% filter(b_scaled_herbivory > 0) %>% nrow /
mod_cwmn %>% spread_draws(b_scaled_herbivory) %>% nrow
## [1] 0.99825
ce_mod_cwmn <- conditional_effects(mod_cwmn, "d_herbivory")
# regression
plot(
ce_mod_cwmn,
plot = FALSE,
points = TRUE,
point_args = list(
fill = "#00ff40",
alpha = 0.4,
pch = 21,
size = 3
),
line_args = list(col = "black", lwd = .6, lty = 1)
)[[1]] +
scale_fill_manual(values = c("#00ff40")) +
scale_colour_manual(values = c("#00ff40")) +
geom_ribbon(fill = "#00ff40", alpha = 0.3) +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("\u0394 CWM N-number") +
xlab("") -> fig4a1
# boxplot
ggplot(data = d3,
aes(y = d_cwmn, x = 1)) +
geom_violinhalf(fill = "#00ff40", alpha = 0.4) +
geom_point(
aes(y = d_cwmn, x = 0.92),
position = position_jitter(width = 0.04),
col = "#00ff40",
alpha = 0.3
) +
geom_boxplot(
aes(x = 0.92) ,
outlier.shape = NA,
width = 0.1,
alpha = 0.4,
fill = "#00ff40"
) +
stat_summary(
fun = mean,
geom = "point",
shape = 24,
size = 4,
color = "black",
fill = "#00ff40",
alpha = 0.8
) +
geom_hline(yintercept = 0) +
labs(x = "", y = "") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.text.x = element_blank(),
axis.text.y = element_blank(),
aspect.ratio = 2
) +
scale_y_continuous(position = "right") -> fig4a2
fig4a1 + fig4a2
mod_alien <- brms::brm(data = d3,
scale(d_alien_percent) ~
scale(d_herbivory) +
scale(baseline_herbivory) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 2000, chains = 4, cores = 4,
file = "models/mod_alien")
mod_alien
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(d_alien_percent) ~ scale(d_herbivory) + scale(baseline_herbivory) + scale(time_span) + scale(site_area_log)
## Data: d3 (Number of observations: 52)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -0.00 0.14 -0.27 0.27 1.00 6213
## scaled_herbivory 0.36 0.16 0.04 0.68 1.00 4015
## scalebaseline_herbivory -0.09 0.14 -0.36 0.18 1.00 4854
## scaletime_span 0.04 0.15 -0.25 0.34 1.00 4965
## scalesite_area_log 0.02 0.16 -0.30 0.32 1.00 4943
## Tail_ESS
## Intercept 2898
## scaled_herbivory 2899
## scalebaseline_herbivory 3048
## scaletime_span 3000
## scalesite_area_log 3384
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.97 0.10 0.79 1.19 1.00 5045 3015
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for positive slope
mod_alien %>% spread_draws(b_scaled_herbivory) %>% filter(b_scaled_herbivory > 0) %>% nrow /
mod_alien %>% spread_draws(b_scaled_herbivory) %>% nrow
## [1] 0.98325
ce_mod_alien <- conditional_effects(mod_alien, "d_herbivory")
# regression
plot(
ce_mod_alien,
plot = FALSE,
points = TRUE,
point_args = list(
fill = "#00bfff",
alpha = 0.4,
pch = 21,
size = 3
),
line_args = list(col = "black", lwd = .6, lty = 1)
)[[1]] +
scale_fill_manual(values = c("#00bfff")) +
scale_colour_manual(values = c("#00bfff")) +
geom_ribbon(fill = "#00bfff", alpha = 0.3) +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("\u0394 % non-native spp.") +
xlab("") -> fig4b1
# boxplot
ggplot(data = d3,
aes(y = d_alien_percent * 100, x = 1)) +
geom_violinhalf(fill = "#00bfff", alpha = 0.4) +
geom_point(
aes(y = d_alien_percent * 100, x = 0.92),
position = position_jitter(width = 0.04),
col = "#00bfff",
alpha = 0.3
) +
geom_boxplot(
aes(x = 0.92) ,
outlier.shape = NA,
width = 0.1,
alpha = 0.4,
fill = "#00bfff"
) +
stat_summary(
fun = mean,
geom = "point",
shape = 24,
size = 4,
color = "black",
fill = "#00bfff",
alpha = 0.8
) +
geom_hline(yintercept = 0) +
labs(x = "", y = "") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.text.x = element_blank(),
axis.text.y = element_blank(),
aspect.ratio = 2
) +
scale_y_continuous(position = "right") -> fig4b2
fig4b1 + fig4b2
mod_redlist <- brms::brm(data = d3,
scale(d_redlist_percent) ~
scale(d_herbivory) +
scale(baseline_herbivory) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 2000, chains = 4, cores = 4,
file = "models/mod_redlist")
mod_redlist
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(d_redlist_percent) ~ scale(d_herbivory) + scale(baseline_herbivory) + scale(time_span) + scale(site_area_log)
## Data: d3 (Number of observations: 52)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 0.00 0.14 -0.27 0.27 1.00 5316
## scaled_herbivory -0.19 0.17 -0.53 0.14 1.00 3741
## scalebaseline_herbivory -0.06 0.15 -0.37 0.23 1.00 5652
## scaletime_span -0.06 0.16 -0.37 0.25 1.00 4619
## scalesite_area_log -0.14 0.16 -0.46 0.18 1.00 4807
## Tail_ESS
## Intercept 3234
## scaled_herbivory 2895
## scalebaseline_herbivory 3068
## scaletime_span 3114
## scalesite_area_log 3179
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.01 0.11 0.82 1.25 1.00 4577 2872
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for negative slope
mod_redlist %>% spread_draws(b_scaled_herbivory) %>% filter(b_scaled_herbivory < 0) %>% nrow /
mod_redlist %>% spread_draws(b_scaled_herbivory) %>% nrow
## [1] 0.881
ce_mod_redlist <- conditional_effects(mod_redlist, "d_herbivory")
# regression
plot(
ce_mod_redlist,
plot = FALSE,
points = TRUE,
point_args = list(
fill = "#ff00bf",
alpha = 0.4,
pch = 21,
size = 3
),
line_args = list(col = "black", lwd = .6, lty = 2)
)[[1]] +
scale_fill_manual(values = c("#ff00bf")) +
scale_colour_manual(values = c("##ff00bf")) +
geom_ribbon(fill = "#ff00bf", alpha = 0.3) +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("\u0394 % red-listed spp.") +
xlab("") -> fig4c1
# boxplot
ggplot(data = d3,
aes(y = d_redlist_percent * 100, x = 1)) +
geom_violinhalf(fill = "#ff00bf", alpha = 0.4) +
geom_point(
aes(y = d_redlist_percent * 100, x = 0.92),
position = position_jitter(width = 0.04),
col = "#ff00bf",
alpha = 0.3
) +
geom_boxplot(
aes(x = 0.92) ,
outlier.shape = NA,
width = 0.1,
alpha = 0.4,
fill = "#ff00bf"
) +
stat_summary(
fun = mean,
geom = "point",
shape = 24,
size = 4,
color = "black",
fill = "#ff00bf",
alpha = 0.8
) +
geom_hline(yintercept = 0) +
labs(x = "", y = "") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.text.x = element_blank(),
axis.text.y = element_blank(),
aspect.ratio = 2
) +
scale_y_continuous(position = "right") -> fig4c2
fig4c1 + fig4c2
mod_smallr <- brms::brm(data = d3,
scale(d_small_ranged_percent) ~
scale(d_herbivory) +
scale(baseline_herbivory) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 2000, chains = 4, cores = 4,
file = "models/mod_smallr")
mod_smallr
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(d_small_ranged_percent) ~ scale(d_herbivory) + scale(baseline_herbivory) + scale(time_span) + scale(site_area_log)
## Data: d3 (Number of observations: 52)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -0.00 0.14 -0.28 0.27 1.00 5560
## scaled_herbivory -0.24 0.17 -0.58 0.08 1.00 4456
## scalebaseline_herbivory -0.18 0.15 -0.47 0.11 1.00 4360
## scaletime_span -0.05 0.16 -0.36 0.25 1.00 4832
## scalesite_area_log -0.13 0.17 -0.45 0.20 1.00 4541
## Tail_ESS
## Intercept 2942
## scaled_herbivory 3383
## scalebaseline_herbivory 2927
## scaletime_span 3123
## scalesite_area_log 3008
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.00 0.10 0.82 1.22 1.00 4673 3335
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for negative slope
mod_smallr %>% spread_draws(b_scaled_herbivory) %>% filter(b_scaled_herbivory < 0) %>% nrow /
mod_smallr %>% spread_draws(b_scaled_herbivory) %>% nrow
## [1] 0.92325
ce_mod_smallr <- conditional_effects(mod_smallr, "d_herbivory")
# regression
plot(
ce_mod_smallr,
plot = FALSE,
points = TRUE,
point_args = list(
fill = "#ff4000",
alpha = 0.4,
pch = 21,
size = 3
),
line_args = list(col = "black", lwd = .6, lty = 2)
)[[1]] +
scale_fill_manual(values = c("#ff4000")) +
scale_colour_manual(values = c("#ff4000")) +
geom_ribbon(fill = "#ff4000", alpha = 0.3) +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("\u0394 % small-ranged spp.") +
xlab("\u0394 herbivory pressure") -> fig4d1
# boxplot
ggplot(data = d3,
aes(y = d_small_ranged_percent * 100, x = 1)) +
geom_violinhalf(fill = "#ff4000", alpha = 0.4) +
geom_point(
aes(y = d_small_ranged_percent * 100, x = 0.92),
position = position_jitter(width = 0.04),
col = "#ff4000",
alpha = 0.3
) +
geom_boxplot(
aes(x = 0.92) ,
outlier.shape = NA,
width = 0.1,
alpha = 0.4,
fill = "#ff4000"
) +
stat_summary(
fun = mean,
geom = "point",
shape = 24,
size = 4,
color = "black",
fill = "#ff4000",
alpha = 0.8
) +
geom_hline(yintercept = 0) +
labs(x = "Density", y = "") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.text.x = element_text(colour = 'white'),
axis.text.y = element_text(colour = 'white'),
aspect.ratio = 2
) +
scale_y_continuous(position = "right") -> fig4d2
fig4d1 + fig4d2
mod_cwmn_inter <- brms::brm(data = d3,
scale(d_cwmn) ~
scale(d_herbivory) * scale(cum_ndepo) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 4000, chains = 4, cores = 4,
file = "models/mod_cwmn_inter")
mod_cwmn_inter
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(d_cwmn) ~ scale(d_herbivory) * scale(cum_ndepo) + scale(time_span) + scale(site_area_log)
## Data: d3 (Number of observations: 52)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept -0.14 0.14 -0.41 0.12 1.00
## scaled_herbivory 0.28 0.16 -0.02 0.59 1.00
## scalecum_ndepo 0.24 0.20 -0.14 0.63 1.00
## scaletime_span 0.10 0.17 -0.23 0.44 1.00
## scalesite_area_log 0.12 0.14 -0.16 0.40 1.00
## scaled_herbivory:scalecum_ndepo 0.33 0.17 0.01 0.66 1.00
## Bulk_ESS Tail_ESS
## Intercept 8097 5979
## scaled_herbivory 6493 6007
## scalecum_ndepo 5215 5302
## scaletime_span 6318 6002
## scalesite_area_log 6944 6112
## scaled_herbivory:scalecum_ndepo 6086 5299
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.84 0.09 0.69 1.04 1.00 7372 5918
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for positive interaction
mod_cwmn_inter %>%
spread_draws(`b_scaled_herbivory:scalecum_ndepo`) %>%
filter(`b_scaled_herbivory:scalecum_ndepo` > 0) %>% nrow /
mod_cwmn_inter %>% spread_draws(`b_scaled_herbivory:scalecum_ndepo`) %>% nrow
## [1] 0.976625
# 10th percentile N deposition
quantile(d3$cum_ndepo, probs = c(0.1, 0.9))
## 10% 90%
## 347.985 1010.374
# counterfactual for low N deposition
int_conditions10 <-
list(cum_ndepo = setNames(c(348), c("low (10th percentile kg/ha)")))
ce_mod_cwmn_inter <- conditional_effects(mod_cwmn_inter,
"d_herbivory:cum_ndepo",
int_conditions = int_conditions10)
plot(
ce_mod_cwmn_inter,
plot = FALSE,
points = FALSE,
point_args = list(col = "grey30", alpha = 0.2),
line_args = list(lwd = .6, alpha = 0.2),
rug = TRUE,
rug_args = list(alpha = 0.8, col = "magenta")
)[[1]] +
scale_fill_manual(values = c("#00ff40"), name = "Cumulative N-dep") +
scale_colour_manual(values = "black") +
labs(y = "\u0394 CWM N-number", x = "", title = "Low N-deposition") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica",
plot_title_face = "bold"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
plot.title = element_text(hjust = 0.5),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) -> fig4a3
# counterfactual for high N deposition
int_conditions90 <-
list(cum_ndepo = setNames(c(1010), c("high (90th percentile kg/ha)")))
ce_mod_cwmn_inter <- conditional_effects(mod_cwmn_inter,
"d_herbivory:cum_ndepo",
int_conditions = int_conditions90)
plot(
ce_mod_cwmn_inter,
plot = FALSE,
points = FALSE,
point_args = list(col = "grey30", alpha = 0.6),
line_args = list(lwd = .6, alpha = 0.6),
rug = TRUE,
rug_args = list(alpha = 0.8, col = "magenta")
)[[1]] +
scale_fill_manual(values = c("#00ff40"), name = "Cumulative N-dep") +
scale_colour_manual(values = "black") +
labs(y = "", x = "", title = "High N-deposition") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica",
plot_title_face = "bold"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
plot.title = element_text(hjust = 0.5),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) -> fig4a4
fig4a3 + fig4a4
mod_alien_inter <- brms::brm(data = d3,
scale(d_alien_percent) ~
scale(d_herbivory) * scale(cum_ndepo) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 4000, chains = 4, cores = 4,
file = "models/mod_alien_inter")
mod_alien_inter
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(d_alien_percent) ~ scale(d_herbivory) * scale(cum_ndepo) + scale(time_span) + scale(site_area_log)
## Data: d3 (Number of observations: 52)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept -0.16 0.14 -0.44 0.13 1.00
## scaled_herbivory 0.19 0.16 -0.12 0.51 1.00
## scalecum_ndepo 0.66 0.21 0.26 1.07 1.00
## scaletime_span -0.26 0.18 -0.62 0.08 1.00
## scalesite_area_log -0.07 0.15 -0.38 0.23 1.00
## scaled_herbivory:scalecum_ndepo 0.36 0.17 0.01 0.70 1.00
## Bulk_ESS Tail_ESS
## Intercept 7862 5991
## scaled_herbivory 6429 6270
## scalecum_ndepo 5601 5531
## scaletime_span 6151 5508
## scalesite_area_log 6498 5817
## scaled_herbivory:scalecum_ndepo 6748 6149
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.88 0.10 0.71 1.09 1.00 6816 5916
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for positive interaction
mod_alien_inter %>%
spread_draws(`b_scaled_herbivory:scalecum_ndepo`) %>%
filter(`b_scaled_herbivory:scalecum_ndepo` > 0) %>% nrow /
mod_alien_inter %>% spread_draws(`b_scaled_herbivory:scalecum_ndepo`) %>% nrow
## [1] 0.97925
# counterfactual for low N deposition
ce_mod_alien_inter <- conditional_effects(mod_alien_inter,
"d_herbivory:cum_ndepo",
int_conditions = int_conditions10)
plot(
ce_mod_alien_inter,
plot = FALSE,
points = FALSE,
point_args = list(col = "grey30", alpha = 0.25),
line_args = list(lwd = .6, alpha = 0.25),
rug = TRUE,
rug_args = list(alpha = 0.8, col = "magenta")
)[[1]] +
scale_fill_manual(values = c("#00bfff"), name = "Cumulative N-dep") +
scale_colour_manual(values = "black") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("\u0394 % non-native spp.") +
xlab("") -> fig4b3
# counterfactual for high N deposition
ce_mod_alien_inter <- conditional_effects(mod_alien_inter,
"d_herbivory:cum_ndepo",
int_conditions = int_conditions90)
plot(
ce_mod_alien_inter,
plot = FALSE,
points = FALSE,
point_args = list(col = "grey30", alpha = 0.6),
line_args = list(lwd = .6, alpha = 0.5),
rug = TRUE,
rug_args = list(alpha = 0.8, col = "magenta")
)[[1]] +
scale_fill_manual(values = c("#00bfff"), name = "Cumulative N-dep") +
scale_colour_manual(values = "black") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("") +
xlab("") -> fig4b4
fig4b3 + fig4b4
mod_redlist_inter <- brms::brm(data = d3,
scale(d_redlist_percent) ~
scale(d_herbivory) * scale(cum_ndepo) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 4000, chains = 4, cores = 4,
file = "models/mod_redlist_inter")
mod_redlist_inter
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(d_redlist_percent) ~ scale(d_herbivory) * scale(cum_ndepo) + scale(time_span) + scale(site_area_log)
## Data: d3 (Number of observations: 52)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 0.16 0.15 -0.14 0.47 1.00
## scaled_herbivory 0.01 0.18 -0.33 0.36 1.00
## scalecum_ndepo -0.60 0.22 -1.05 -0.17 1.00
## scaletime_span 0.19 0.19 -0.17 0.57 1.00
## scalesite_area_log -0.06 0.16 -0.39 0.26 1.00
## scaled_herbivory:scalecum_ndepo -0.38 0.19 -0.75 -0.01 1.00
## Bulk_ESS Tail_ESS
## Intercept 8236 6186
## scaled_herbivory 5880 5837
## scalecum_ndepo 5094 4814
## scaletime_span 6349 5910
## scalesite_area_log 6459 5724
## scaled_herbivory:scalecum_ndepo 6144 5567
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.94 0.10 0.77 1.15 1.00 6749 5773
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for negative interaction
mod_redlist_inter %>%
spread_draws(`b_scaled_herbivory:scalecum_ndepo`) %>%
filter(`b_scaled_herbivory:scalecum_ndepo` < 0) %>% nrow /
mod_redlist_inter %>% spread_draws(`b_scaled_herbivory:scalecum_ndepo`) %>% nrow
## [1] 0.978875
# counterfactual for low N deposition
ce_mod_redlist_inter <- conditional_effects(mod_redlist_inter,
"d_herbivory:cum_ndepo",
int_conditions = int_conditions10)
plot(
ce_mod_redlist_inter,
plot = FALSE,
points = FALSE,
point_args = list(col = "grey30", alpha = 0.3),
line_args = list(lwd = .6, alpha = 0.25),
rug = TRUE,
rug_args = list(alpha = 0.8, col = "magenta")
)[[1]] +
scale_fill_manual(values = c("#ff00bf"), name = "Cumulative N-dep") +
scale_colour_manual(values = c("black")) +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("\u0394 % red-listed spp.") +
xlab("") -> fig4c3
# counterfactual for high N deposition
ce_mod_redlist_inter <- conditional_effects(mod_redlist_inter,
"d_herbivory:cum_ndepo",
int_conditions = int_conditions90)
plot(
ce_mod_redlist_inter,
plot = FALSE,
points = FALSE,
point_args = list(col = "grey30", alpha = 0.6),
line_args = list(lwd = .6, alpha = 0.5),
rug = TRUE,
rug_args = list(alpha = 0.8, col = "magenta")
)[[1]] +
scale_fill_manual(values = c("#ff00bf"), name = "Cumulative N-dep") +
scale_colour_manual(values = c("black")) +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("") +
xlab("") -> fig4c4
fig4c3 + fig4c4
mod_smallr_inter <- brms::brm(data = d3,
scale(d_small_ranged_percent) ~
scale(d_herbivory) * scale(cum_ndepo) +
scale(time_span) +
scale(site_area_log),
family = gaussian, iter = 4000, chains = 4, cores = 4,
file = "models/mod_smallr_inter")
mod_smallr_inter
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: scale(d_small_ranged_percent) ~ scale(d_herbivory) * scale(cum_ndepo) + scale(time_span) + scale(site_area_log)
## Data: d3 (Number of observations: 52)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 0.22 0.15 -0.08 0.52 1.00
## scaled_herbivory 0.04 0.17 -0.30 0.38 1.00
## scalecum_ndepo -0.57 0.21 -0.99 -0.15 1.00
## scaletime_span 0.15 0.19 -0.22 0.53 1.00
## scalesite_area_log -0.08 0.16 -0.39 0.23 1.00
## scaled_herbivory:scalecum_ndepo -0.51 0.18 -0.87 -0.15 1.00
## Bulk_ESS Tail_ESS
## Intercept 7357 5798
## scaled_herbivory 6282 5431
## scalecum_ndepo 5339 4777
## scaletime_span 6175 5215
## scalesite_area_log 6420 5019
## scaled_herbivory:scalecum_ndepo 6829 5876
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.92 0.10 0.75 1.14 1.00 5687 5522
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# calculate posterior probability mass for negative interaction
mod_smallr_inter %>%
spread_draws(`b_scaled_herbivory:scalecum_ndepo`) %>%
filter(`b_scaled_herbivory:scalecum_ndepo` < 0) %>% nrow /
mod_smallr_inter %>% spread_draws(`b_scaled_herbivory:scalecum_ndepo`) %>% nrow
## [1] 0.9965
# counterfactual for low N deposition
ce_mod_smallr_inter <- conditional_effects(mod_smallr_inter,
"d_herbivory:cum_ndepo",
int_conditions = int_conditions10)
plot(
ce_mod_smallr_inter,
plot = FALSE,
points = FALSE,
point_args = list(col = "grey30", alpha = 0.3),
line_args = list(lwd = .6, alpha = 0.25),
rug = TRUE,
rug_args = list(alpha = 0.8, col = "magenta")
)[[1]] +
scale_fill_manual(values = c("#ff4000"), name = "Cumulative N-dep") +
scale_colour_manual(values = "black") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("\u0394 % small-ranged spp.") +
xlab("\u0394 herbivory pressure") -> fig4d3
# counterfactual for high N deposition
ce_mod_smallr_inter <- conditional_effects(mod_smallr_inter,
"d_herbivory:cum_ndepo",
int_conditions = int_conditions90)
plot(
ce_mod_smallr_inter,
plot = FALSE,
points = FALSE,
point_args = list(col = "grey30", alpha = 0.65),
line_args = list(lwd = .6, alpha = 0.5),
rug = TRUE,
rug_args = list(alpha = 0.8, col = "magenta")
)[[1]] +
scale_fill_manual(values = c("#ff4000"), name = "Cumulative N-dep") +
scale_colour_manual(values = "black") +
theme_ipsum(
grid = "",
axis_title_size = 12,
plot_title_size = 12,
axis_title_just = "mm",
base_family = "Helvetica"
) +
theme(
plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.position = "none",
axis.line = element_line(),
aspect.ratio = 1
) +
ylab("") +
xlab("\u0394 herbivory pressure") -> fig4d4
fig4d3 + fig4d4
fig4a1 + fig4a2 + fig4a3 + fig4a4 +
fig4b1 + fig4b2 + fig4b3 + fig4b4 +
fig4c1 + fig4c2 + fig4c3 + fig4c4 +
fig4d1 + fig4d2 + fig4d3 + fig4d4 +
plot_layout(ncol = 4,
widths = c(1, 0.5, 1, 1)) +
plot_annotation(tag_levels = list(c("a", "", "b", "",
"c", "", "d", "",
"e", "", "f", "",
"g", "", "h", ""))) &
theme(plot.tag = element_text(face = 'bold'))