Supplementary note: R code for statistical analyses
Josiane Segar and Ingmar Staude
- Preamble
- Load data
- Analyses and visualization
- Shrub, herb and tree layer changes versus changes in herbivory
- Species richness change and exchange ratio versus changes in herbivory
- Changes in community composition versus changes in herbivory
- Model: Community weighted mean (CWM) N-number
- Plot: Community weighted mean (CWM) N-number
- Model: Percent change in non-native species
- Plot: Percent change in non-native species
- Model: Percent change in species listed on national Red Lists
- Plot: Percent change in species listed on national Red Lists
- Model: Percent change in small-ranged species
- Plot: Percent change in small-ranged species
- Changes in community composition versus changes in herbivory in relation to N deposition
- Model: Community weighted mean (CWM) N-number, interaction
- Plot: Community weighted mean (CWM) N-number, interaction
- Model: Percent change in non-native species, interaction
- Plot: Percent change in non-native species, interaction
- Model: Percent change in species listed on national Red Lists, interaction
- Plot: Percent change in species listed on national Red Lists, interaction
- Model: Percent change in small-ranged species, interaction
- Plot: Percent change in small-ranged species, interaction
- Multipanel Figure 4
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:
- id = the site ids from the forestREplot database
- d_herbivory = the change in herbivore pressure over time
- baseline_herbivory = baseline herbivore pressure
- time_span = the time interval between baseline and resurvey
- site_area_log = the total site area of a study site
- herb_cover = the change in herb cover over time
- shrub_cover = the change in shrub cover over time
- tree_cover = the change in tree cover over time
- d_spp_richness = the change in species richness over time
- spp_exchange_ratio = species richness exchange ratio
- d_small_ranged_percent = the percentage change in small-ranged species over time
- d_redlist_percent = the percentage change in red-listed species over time
- d_cwmn = the change in the CWM N-number of a site over time
- d_alien_percent = the percentage change in non-native species over time
- cum_ndepo = cumulative N deposition on a site between baseline and resurvey
Preamble
# load packages
library(tidyverse)
library(tidybayes)
library(brms)
library(patchwork)
library(hrbrthemes)
library(see)
Load data
# 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
Analyses and visualization
Shrub, herb and tree layer changes versus changes in herbivory
Model: shrub cover change
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
Plot: shrub cover change
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") -> fig2a2fig2a1 + fig2a2
Model: herb cover change
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
Plot: herb cover change
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") -> fig2b2fig2b1 + fig2b2
Model: tree cover change
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
Plot: tree cover change
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") -> fig2c2fig2c1 + fig2c2
Multipanel Figure 2
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')) 
Species richness change and exchange ratio versus changes in herbivory
Model: species richness
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
Plot: species richness
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") -> fig3a2fig3a1 + fig3a2
Model: species exchange ratio
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
Plot: species exchange ratio
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") -> fig3b2fig3b1 + fig3b2
Multipanel Figure 3
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')) 
Changes in community composition versus changes in herbivory
Model: Community weighted mean (CWM) N-number
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
Plot: Community weighted mean (CWM) N-number
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") -> fig4a2fig4a1 + fig4a2
Model: Percent change in non-native species
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
Plot: Percent change in non-native species
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") -> fig4b2fig4b1 + fig4b2
Model: Percent change in species listed on national Red Lists
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
Plot: Percent change in species listed on national Red Lists
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") -> fig4c2fig4c1 + fig4c2
Model: Percent change in small-ranged species
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
Plot: Percent change in small-ranged species
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") -> fig4d2fig4d1 + fig4d2
Changes in community composition versus changes in herbivory in relation to N deposition
Model: Community weighted mean (CWM) N-number, interaction
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
Plot: Community weighted mean (CWM) N-number, interaction
# 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
) -> fig4a4fig4a3 + fig4a4
Model: Percent change in non-native species, interaction
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
Plot: Percent change in non-native species, interaction
# 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("") -> fig4b4fig4b3 + fig4b4
Model: Percent change in species listed on national Red Lists, interaction
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
Plot: Percent change in species listed on national Red Lists, interaction
# 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("") -> fig4c4fig4c3 + fig4c4
Model: Percent change in small-ranged species, interaction
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
Plot: Percent change in small-ranged species, interaction
# 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") -> fig4d4fig4d3 + fig4d4
Multipanel Figure 4
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')) 