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:

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") -> fig2a2
fig2a1 + 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") -> fig2b2
fig2b1 + 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") -> fig2c2
fig2c1 + 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") -> fig3a2
fig3a1 + 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") -> fig3b2
fig3b1 + 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") -> fig4a2
fig4a1 + 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") -> fig4b2
fig4b1 + 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") -> fig4c2
fig4c1 + 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") -> fig4d2
fig4d1 + 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
  ) -> fig4a4
fig4a3 + 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("") -> fig4b4
fig4b3 + 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("") -> fig4c4
fig4c3 + 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") -> fig4d4
fig4d3 + 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'))