multiplestressR - R Package
Introductory Tutorial
Introduction to multiplestressR
One of the common aims of nearly all studies within multiple stressor ecology is to determine how co-occurring stressors combine to affect individuals, populations, or even entire ecological communities. Often, it is simply assumed that where more than one stressor affects a species, the combined effect of these stressors upon the species will simply be the sum of the individual effects. For instance, in the absence of any stressors there may be 100 individuals of Daphnia pulex. However, in the presence of a pesticide there may only 75 individuals, while an increase in temperature may reduce to population size to only 60 individuals. As such, if we sum the effects of the individual stressors (i.e., 100 + (75 - 100) + (60 - 100)), we would assume that in the presence of the pesticide and increased temperature, the population of Daphnia pulex would be reduced to 35 individuals. This expectation of stressor effects (based on an assumption of the additivity of their effects) is refered to as the additive null model. If this assumption is found to be met, then these stressors may be found as having an additive or null interaction.
However, it may be that in the presence of both stressors this assumption is not met. If the population size of Daphnia pulex is found to have been reduced by more than expected than the additive null model (e.g., only 20 individuals survive) then a synergistic interaction is said to be occurring. If the population size is reduced by less than expected by the additive null model (e.g., 50 individuals survive) then an antagonistic interaction is said to be occurring. In some cases it may be that stressors have an effect which diametrically opposes than expected by the additive null model (e.g., the null model predicts a decline in a population size, but under both stressors the population size actually increases). In such a situation a reversal interaction is said to be occurring. However, sometimes no distinction is made between antagonisitc and reversal interactions, with both being described as being an antagonistic interaction.
Understanding how multiple stressors interact to affect ecosystems is crucial. Indeed, knowledge of stressor interactions is crucial when implementing conservation measures, as a failure to consider how stressors interact may render these actions ineffective. However, understading how stressors interact is very rarely easy. Ecological data is frequently noisy (i.e., high levels of variation) while ecological experiments often times have low numbers of replicates, meaning that the statistical power of any experiment may be low. As such, given these factors, the simple summing approach used in the above example is inappropriate for ecological data. Instead more complex statistical methods are required. These methods often require experiment to be conducted with four treatments in a factorial design (i.e., i) control conditions, ii) only stressor A, iii) only stressor B, and iv) both stressors A and B). According for each of these treatments, means (e.g., population densities, or survival rates), uncertainty (e.g., standard deviations), and the number of replicates per treatment are required.
With data on means, uncertainty, and numbers of replciates for each of the four treatments, it is possible to classify stressor interactions using either of the following approaches. The additive null model is implemented through the use of Hedges’ d, while the multiplicative null model is implemented through the use of the response ratio. These statistical tools have been thoroughly documented elsewhere, and so the equations underpinning these null models will not be discussed more within this tutorial. However,for more information on these statistical tools see Burgess et al. (2021a).
However, while these detailed and rigorous methods allow researchers to classify interactions, many previous studies have implemented different or untested (and in some cases incorrect) versions of these null models. As such, these differences in methodologies may be responsible for the inconsistent results across studies in multiple stressor ecology.
Accordingly, the multiplestressR R package has been developed to provide a simple way for researchers to implement either the additive or multiplicative null models upon their datasets. The package implements statistically rigorous versions of these null models, and can be conducted in a few lines of code. As such, this package allows those with minimal knowledge of R, or a limited understanding of the null models to statistically analyse their data. The multiplestressR package is appropriate for those researchers analysing data for single experiments, or those conducting meta-analyses.
Outlined below is a step-by-step tutorial for analysing multiple stressor data using the multiplestressR package.
Throughout this tutorial, R code is outlined in blue, successful output in green, and errors (or similar) in red.
Installation
The first thing that needs to be done is to install the multiplestressR package. This can be installed from CRAN using the following code:
#install.packages("multiplestressR")
##Note that to install the multiplestressR package from CRAN, the above line must be uncommented.
library(multiplestressR)Barring any issues, the package is now installed.
Data
The next step is to load data into R. For this tutorial, we will be using an example dataset which is provided by the multiplestressR package, which can be loaded using the following code:
df <- multiplestressR::survivalWe can view the first few rows of this dataframe using the following code:
head(df)## Sample_Size_Control Standard_Deviation_Control Mean_Control
## 1 4 0.02886751 0.9550000
## 2 10 0.08743887 0.8430000
## 3 4 0.02645751 0.9550000
## 4 6 0.04195235 0.8900000
## 5 4 0.09055385 0.8600000
## 6 6 0.04676181 0.9033333
## Sample_Size_Temperature Standard_Deviation_Temperature Mean_Temperature
## 1 4 0.07023769 0.7700000
## 2 10 0.11246728 0.8560000
## 3 4 0.10424331 0.7800000
## 4 6 0.07174027 0.8066667
## 5 4 0.09814955 0.7750000
## 6 6 0.11338724 0.7783333
## Sample_Size_pH Standard_Deviation_pH Mean_pH Sample_Size_Temperature_pH
## 1 4 0.04041452 0.8950000 4
## 2 10 0.08753412 0.7880000 10
## 3 4 0.03201562 0.7675000 4
## 4 6 0.04578937 0.8683333 6
## 5 4 0.04966555 0.7300000 4
## 6 6 0.04578937 0.9016667 6
## Standard_Deviation_Temperature_pH Mean_Temperature_pH
## 1 0.08770215 0.6725
## 2 0.07397447 0.7750
## 3 0.07762087 0.6575
## 4 0.07563068 0.6900
## 5 0.07365460 0.7275
## 6 0.04037326 0.8450This dataset has been generated (i.e., does not come from ‘real-world’ experiments), but reflects survival data for 250 unique experiments. For this dataset, survival rates are measured for organisms exposed to i) control conditions, ii) altered temperature, iii) altered pH, or iv) altered temperature and pH together.
As you can see, the dataframe contains data for eight interactions (rows) and twelve variables (columns), with three columns for each treatment (X_Control, X_Temperature, X_pH, X_Temperature_pH). For each treatment, there is a column for treatment means (Mean_X), treatment standard deviation (Standard_Deviation_X), and treatment sample size (Sample_Size_X). It is important to note that any dataset must contain these three variables; means must be used (i.e., not medians or other metrics), standard deviation must be reported (i.e., not standard error or confidence intervals), and sample sizes must be reported.
Null models
The next step is to calculate the null models for each interaction.
Within the multiplestressR package, there are two functions for calculating the null models, either effect_size_additive or effect_size_multiplicative. The former is the implementing the additive null model, while the latter is for implementing the multiplicative null model. For this tutorial we will be implementing the additive null model. Though an example of code using the multiplicative null model is detailed at the very end of this tutorial.
For the effect_size_additive function, fourteen variables need to be specified. Of these twelve relate to the previous discussed means, standard deviations, and sample sizes. Specifying the variables in this manner ensures that each treatment is assigned the correct data (e.g., Control and StressorA are not mixed up) and that different variables are not incorrectly labelled (e.g., standard deviations and means are not mixed up).
The final two variables are user specified at this point. Firstly, Significance_Level corresponds to the level of alpha which should be used when calculating confidence intervals for a given interaction effect size. The default is 0.05, (i.e., calculating 95% confidence intervals) but this can be changed by specifying a number between 0 and 1. Secondly, Hedges’ d (i.e., the form of the additive null model used here) has been shown to slightly overestimate effects where small sizes are small (i.e., <20). As such, there is a small statistical correction which can be applied to the null model to overcome this bias. By default, the variable Small_Sample_Correction is set to TRUE, namely the bias is corrected, though this can be disabled by setting Small_Sample_Correction to FALSE. Note that for the mulitplicative null model (i.e., effect_size_multiplicative this bias does not exist and so the parameter is specified in the model).
The additive null model can be run using the following code:
df <- effect_size_additive(Control_N = df$Sample_Size_Control,
Control_SD = df$Standard_Deviation_Control,
Control_Mean = df$Mean_Control,
StressorA_N = df$Sample_Size_Temperature,
StressorA_SD = df$Standard_Deviation_Temperature,
StressorA_Mean = df$Mean_Temperature,
StressorB_N = df$Sample_Size_pH,
StressorB_SD = df$Standard_Deviation_pH,
StressorB_Mean = df$Mean_pH,
StressorsAB_N = df$Sample_Size_Temperature_pH,
StressorsAB_SD = df$Standard_Deviation_Temperature_pH,
StressorsAB_Mean = df$Mean_Temperature_pH,
Small_Sample_Correction = TRUE,
Significance_Level = 0.05)
head(df)## Control_N Control_SD Control_Mean StressorA_N StressorA_SD StressorA_Mean
## 1 4 0.02886751 0.9550000 4 0.07023769 0.7700000
## 2 10 0.08743887 0.8430000 10 0.11246728 0.8560000
## 3 4 0.02645751 0.9550000 4 0.10424331 0.7800000
## 4 6 0.04195235 0.8900000 6 0.07174027 0.8066667
## 5 4 0.09055385 0.8600000 4 0.09814955 0.7750000
## 6 6 0.04676181 0.9033333 6 0.11338724 0.7783333
## StressorB_N StressorB_SD StressorB_Mean StressorsAB_N StressorsAB_SD
## 1 4 0.04041452 0.8950000 4 0.08770215
## 2 10 0.08753412 0.7880000 10 0.07397447
## 3 4 0.03201562 0.7675000 4 0.07762087
## 4 6 0.04578937 0.8683333 6 0.07563068
## 5 4 0.04966555 0.7300000 4 0.07365460
## 6 6 0.04578937 0.9016667 6 0.04037326
## StressorsAB_Mean Interaction_Effect_Size Interaction_Variance
## 1 0.6725 -0.4949691 0.8831246
## 2 0.7750 -0.2641530 0.3842288
## 3 0.6575 0.7724633 0.8927570
## 4 0.6900 -1.3751384 0.6534558
## 5 0.7275 0.8340493 0.8954668
## 6 0.8450 0.8760449 0.6317925
## Interaction_CI_Upper Interaction_CI_Lower Null_Model
## 1 1.3469013 -2.3368395 Additive
## 2 0.9507540 -1.4790601 Additive
## 3 2.6243513 -1.0794248 Additive
## 4 0.2092301 -2.9595069 Additive
## 5 2.6887457 -1.0206471 Additive
## 6 2.4339296 -0.6818399 AdditiveAccordingly, this function adds five new variables (columns) and slightly alters the dataset.
Firstly, note that within the dataset, the columns for means, standard deviations, and sample sizes have been renamed.
Secondly, Interaction_Effect_Size is a value which indicates the strength of any interaction which is occuring, the greater the value (positive or negative) of the Interaction_Effect_Size the stronger the interaction. If an interaction has an Interaction_Effect_Size of exactly 0, then the effect predicted by the additive null model is exactly the same as the observed effect of the interacting stressors.
Thirdly, Interaction_Variance details the uncertainty in the Interaction_Effect_Size. Given the uncertainty (i.e., standard deviations) in the raw data it stands to reason that there will also be uncertainty in the value of Interaction_Effect_Size. This value for the variance can be be used to calculate confidence intervals based upon the chosen Significance_Level here referred to as Interaction_CI_Upper and Interaction_CI_Lower. The values of Interaction_Effect_Size and confidence intervals can then be used to classify each interaction.
The final variable in the dataset is Null_Model, this simply reflects whether the additive or multiplicative null model was implemented.
Classifying interactions
The next step is to classify each interaction.
Within the multiplestressR package, this can be done using the classify_interactions function.
To implement this function, the following three variables need to be specified.
Firstly, effect_size_dataframe is the output of the either the effect_size_additive or effect_size_multiplicative functions. Note that this must be the exact output of either of these functions.
Secondly, assign_reversals determines whether the interaction classificaton of reversal should be distinguished from antagonisms (see above). The default for this function is TRUE; though this can be altered by specifying the variable as FALSE.
Thirdly, the remove_directionality parameter can be implemented. For the purposes of this tutorial this parameter is not considered any further, its default value is likewise FALSE. However, while this parameter is unlikely to be important for the results of any single experiment, it is important for those researchers conducting a meta-analysis. This will be addressed in a subsequent tutorial. At present, those conducting a meta-analysis should consult the following published papers when considering the importance of removing directionality from those analyses (Burgess et al. (2021b), Jackson et al. (2016), Seifert et al. (2020), Siviter et al. (2021)).
For a detailed outline of this function (e.g., how interactions are assigned, and the use of the remove_directionality parameter), see the help guide for the classify_interactions function. This can be accessed by using the following code:
# ?multiplestressR::classify_interactions
## Note that the above `#` must be removed in order for this code to be run.As such, interactions can be classified by using the following code:
df <- classify_interactions(effect_size_dataframe = df,
assign_reversals = TRUE)
head(df)## Control_N Control_SD Control_Mean StressorA_N StressorA_SD StressorA_Mean
## 1 4 0.02886751 0.9550000 4 0.07023769 0.7700000
## 2 10 0.08743887 0.8430000 10 0.11246728 0.8560000
## 3 4 0.02645751 0.9550000 4 0.10424331 0.7800000
## 4 6 0.04195235 0.8900000 6 0.07174027 0.8066667
## 5 4 0.09055385 0.8600000 4 0.09814955 0.7750000
## 6 6 0.04676181 0.9033333 6 0.11338724 0.7783333
## StressorB_N StressorB_SD StressorB_Mean StressorsAB_N StressorsAB_SD
## 1 4 0.04041452 0.8950000 4 0.08770215
## 2 10 0.08753412 0.7880000 10 0.07397447
## 3 4 0.03201562 0.7675000 4 0.07762087
## 4 6 0.04578937 0.8683333 6 0.07563068
## 5 4 0.04966555 0.7300000 4 0.07365460
## 6 6 0.04578937 0.9016667 6 0.04037326
## StressorsAB_Mean Interaction_Effect_Size Interaction_Variance
## 1 0.6725 -0.4949691 0.8831246
## 2 0.7750 -0.2641530 0.3842288
## 3 0.6575 0.7724633 0.8927570
## 4 0.6900 -1.3751384 0.6534558
## 5 0.7275 0.8340493 0.8954668
## 6 0.8450 0.8760449 0.6317925
## Interaction_CI_Upper Interaction_CI_Lower Null_Model
## 1 1.3469013 -2.3368395 Additive
## 2 0.9507540 -1.4790601 Additive
## 3 2.6243513 -1.0794248 Additive
## 4 0.2092301 -2.9595069 Additive
## 5 2.6887457 -1.0206471 Additive
## 6 2.4339296 -0.6818399 Additive
## Interaction_Classification
## 1 Null
## 2 Null
## 3 Null
## 4 Null
## 5 Null
## 6 NullBy considering the first few rows of the dataset, it is evident that a single new column has been added to the dataset. This column Interaction_Classification describes the classification assigned to see interaction.
Summary plots
The multiplestressR package also has a function (summary_plots) which can be used to generate five plots which may help interpret any analyses. This function has three input parameters,
Firstly, effect_size_dataframe which corresponds to the output of the classify_interactions function.
Secondly, Small_Sample_Correction which corresponds to whether or not the bias correction for small sample sizes was implemented in the effect_size_additive function.
Finally, Significance_Level which again corresponds to the significance level used in the effect_size_additive (or effect_size_multiplicative) function.
The summary plots can be generated using the following code:
dfa_plots <- summary_plots(effect_size_dataframe = df,
Small_Sample_Correction = TRUE,
Significance_Level = 0.05)
dfa_plots
For this analysis, only the first two plots (going clockwise from top-left) are overly useful for any interpretation, although the final three plots are likely to be useful to those researchers conducting a meta-analysis.
The first figure reveals the frequency of the different interaction classifications for interactions within a dataset. For this analysis it is apparent that 76.4% of interactions are assigned an null (i.e., additive) interaction class, 16.8% of interactions are assigned a antagonistic class, 4% of interactions are assigned a synergistic interaction class, while 2.8% of interactions are assigned a reversal interaction class.
The second figure plots Interaction_Effect_Size against median treatment sample sizes (averaged across all four treatments). On this figure, black lines indicate the critical effect size for a given sample size (i.e., the minimum effect size required in order for an interaction to be significantly different to zero). Only those interactions which have interaction effect sizes greater than the critical effect size have been assigned a non-null interaction classification.
It is also important to note that a null interaction classification does not necessarily mean that no interaction is occurring. A null classification may instead be assigned due to an experiment having an insufficient power to detect a given interaction (i.e., insufficient sample sizes to statistically determine whether an interaction between stressors was occurring).
For more information on the role of sample sizes in determining interaction classifications and critical effect sizes, see Burgess et al. (2021a).
Finally, all these figures are generated using ggplot, as such they can be altered by the user. For instance, if a user wishes to use the second figure with the legend shown and using a logarithmic scale on the x-axis, this can be done using the following code.
library(ggplot2)
dfa_plots[[2]] +
scale_x_continuous(trans='log10') +
theme(legend.position="top",
legend.title = element_blank())
Concluding remarks
Overall, this tutorial has shown how the multiplestressR package can be easily and quickly used to implement null models and classify interactions. Indeed, it provides an easy to use framework which is applicable in a range of scenarios. Furthermore, this package allows the results of experiments, studies, or meta-analyses to be easily compared with no concerns regarding how differences in analytical methodologies may affect results.
This tutorial likewise provides a framework which researchers may wish to modify for their studies, allowing researchers unfamiliar with R or these statistical tools to easily conduct a rigorous analysis of their data.
Multiplicative null model
Below is an example of the code which can be used to conduct the same analysis as that detailed above, but for the multiplicative null model, not the additive null model.
#install.packages("multiplestressR")
##Note that to install the multiplestressR package from CRAN, the above line must be uncommented.
library(multiplestressR)
dfm <- multiplestressR::survival
head(dfm)## Sample_Size_Control Standard_Deviation_Control Mean_Control
## 1 4 0.02886751 0.9550000
## 2 10 0.08743887 0.8430000
## 3 4 0.02645751 0.9550000
## 4 6 0.04195235 0.8900000
## 5 4 0.09055385 0.8600000
## 6 6 0.04676181 0.9033333
## Sample_Size_Temperature Standard_Deviation_Temperature Mean_Temperature
## 1 4 0.07023769 0.7700000
## 2 10 0.11246728 0.8560000
## 3 4 0.10424331 0.7800000
## 4 6 0.07174027 0.8066667
## 5 4 0.09814955 0.7750000
## 6 6 0.11338724 0.7783333
## Sample_Size_pH Standard_Deviation_pH Mean_pH Sample_Size_Temperature_pH
## 1 4 0.04041452 0.8950000 4
## 2 10 0.08753412 0.7880000 10
## 3 4 0.03201562 0.7675000 4
## 4 6 0.04578937 0.8683333 6
## 5 4 0.04966555 0.7300000 4
## 6 6 0.04578937 0.9016667 6
## Standard_Deviation_Temperature_pH Mean_Temperature_pH
## 1 0.08770215 0.6725
## 2 0.07397447 0.7750
## 3 0.07762087 0.6575
## 4 0.07563068 0.6900
## 5 0.07365460 0.7275
## 6 0.04037326 0.8450dfm <- effect_size_multiplicative(Control_N = dfm$Sample_Size_Control,
Control_SD = dfm$Standard_Deviation_Control,
Control_Mean = dfm$Mean_Control,
StressorA_N = dfm$Sample_Size_Temperature,
StressorA_SD = dfm$Standard_Deviation_Temperature,
StressorA_Mean = dfm$Mean_Temperature,
StressorB_N = dfm$Sample_Size_pH,
StressorB_SD = dfm$Standard_Deviation_pH,
StressorB_Mean = dfm$Mean_pH,
StressorsAB_N = dfm$Sample_Size_Temperature_pH,
StressorsAB_SD = dfm$Standard_Deviation_Temperature_pH,
StressorsAB_Mean = dfm$Mean_Temperature_pH,
Significance_Level = 0.05)
head(dfm)## Control_N Control_SD Control_Mean StressorA_N StressorA_SD StressorA_Mean
## 1 4 0.02886751 0.9550000 4 0.07023769 0.7700000
## 2 10 0.08743887 0.8430000 10 0.11246728 0.8560000
## 3 4 0.02645751 0.9550000 4 0.10424331 0.7800000
## 4 6 0.04195235 0.8900000 6 0.07174027 0.8066667
## 5 4 0.09055385 0.8600000 4 0.09814955 0.7750000
## 6 6 0.04676181 0.9033333 6 0.11338724 0.7783333
## StressorB_N StressorB_SD StressorB_Mean StressorsAB_N StressorsAB_SD
## 1 4 0.04041452 0.8950000 4 0.08770215
## 2 10 0.08753412 0.7880000 10 0.07397447
## 3 4 0.03201562 0.7675000 4 0.07762087
## 4 6 0.04578937 0.8683333 6 0.07563068
## 5 4 0.04966555 0.7300000 4 0.07365460
## 6 6 0.04578937 0.9016667 6 0.04037326
## StressorsAB_Mean Interaction_Effect_Size Interaction_Variance
## 1 0.6725 -0.07050078 0.007070194
## 2 0.7750 -0.03193848 0.004947161
## 3 0.6575 0.04772371 0.008576382
## 4 0.6900 -0.13157314 0.004154372
## 5 0.7275 0.10063882 0.010501229
## 6 0.8450 0.08402847 0.004793997
## Interaction_CI_Upper Interaction_CI_Lower Null_Model
## 1 0.09430170 -0.23530326 Multiplicative
## 2 0.10591767 -0.16979462 Multiplicative
## 3 0.22923354 -0.13378613 Multiplicative
## 4 -0.00524479 -0.25790148 Multiplicative
## 5 0.30148712 -0.10020947 Multiplicative
## 6 0.21973382 -0.05167688 Multiplicative#classifying interactions
dfm <- classify_interactions(effect_size_dataframe = dfm,
assign_reversals = TRUE,
remove_directionality = TRUE)
head(dfm)## Control_N Control_SD Control_Mean StressorA_N StressorA_SD StressorA_Mean
## 1 4 0.02886751 0.9550000 4 0.07023769 0.7700000
## 2 10 0.08743887 0.8430000 10 0.11246728 0.8560000
## 3 4 0.02645751 0.9550000 4 0.10424331 0.7800000
## 4 6 0.04195235 0.8900000 6 0.07174027 0.8066667
## 5 4 0.09055385 0.8600000 4 0.09814955 0.7750000
## 6 6 0.04676181 0.9033333 6 0.11338724 0.7783333
## StressorB_N StressorB_SD StressorB_Mean StressorsAB_N StressorsAB_SD
## 1 4 0.04041452 0.8950000 4 0.08770215
## 2 10 0.08753412 0.7880000 10 0.07397447
## 3 4 0.03201562 0.7675000 4 0.07762087
## 4 6 0.04578937 0.8683333 6 0.07563068
## 5 4 0.04966555 0.7300000 4 0.07365460
## 6 6 0.04578937 0.9016667 6 0.04037326
## StressorsAB_Mean Interaction_Effect_Size Interaction_Variance
## 1 0.6725 0.07050078 0.007070194
## 2 0.7750 0.03193848 0.004947161
## 3 0.6575 -0.04772371 0.008576382
## 4 0.6900 0.13157314 0.004154372
## 5 0.7275 -0.10063882 0.010501229
## 6 0.8450 -0.08402847 0.004793997
## Interaction_CI_Upper Interaction_CI_Lower Null_Model
## 1 0.23530326 -0.09430170 Multiplicative
## 2 0.16979462 -0.10591767 Multiplicative
## 3 0.13378613 -0.22923354 Multiplicative
## 4 0.25790148 0.00524479 Multiplicative
## 5 0.10020947 -0.30148712 Multiplicative
## 6 0.05167688 -0.21973382 Multiplicative
## Interaction_Classification
## 1 Null
## 2 Null
## 3 Null
## 4 Synergistic
## 5 Null
## 6 Null#generate summary plots
dfm_plots <- summary_plots(effect_size_dataframe = dfm,
Significance_Level = 0.05)
dfm_plots
References
Burgess, B. J., Jackson, M. C., & Murrell, D. J. (2021a). Multiple stressor null models frequently fail to detect most interactions due to low statistical power. bioRxiv.
Burgess, B. J., Purves, D., Mace, G., & Murrell, D. J. (2021b). Classifying ecosystem stressor interactions: Theory highlights the data limitations of the additive null model and the difficulty in revealing ecological surprises. Global Change Biology.
Jackson, M. C., Loewen, C. J., Vinebrooke, R. D., & Chimimba, C. T. (2016). Net effects of multiple stressors in freshwater ecosystems: a meta-analysis. Global Change Biology, 22(1), 180-189.
Seifert, M., Rost, B., Trimborn, S., & Hauck, J. (2020). Meta-analysis of multiple driver effects on marine phytoplannkton highlights modulating role of pCO2. Global Change Biology, 26(12), 6787-6804.
Siviter, H., Bailes, E. J., Martin, C. D., Oliver, T. R., Koricheva, J., Leadbeater, E., & Brown, M. J. (2021). Agrochemicals interact synergistically to increase bee mortality. Nature, 596(7872), 389-392.