--- title: "Analysis of scale-dependent biodiversity changes with mobr" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Analysis of scale-dependent biodiversity changes with mobr} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- # Installing `mobr` The package `mobr` is currently available on GitHub and you can freely download the source code here [mobr on GitHub](https://github.com/MoBiodiv/mobr) The easiest option is to install the package directly from GitHub using the package `devtools`. If you do not already have `devtools` installed then need to install it. ```{r install package, eval = F} install.packages('devtools') library(devtools) install_github('MoBiodiv/mobr') ``` If you receive an error we would love to receive your bug report [here](https://github.com/MoBiodiv/mobr/issues) # Data structure required by `mobr` To work with `mobr` you need two matrix-like data tables: 1. Species abundance data in a community matrix (rows are sites and columns are species) 2. A site attributes table (rows are sites and columns are site attributes) The community matrix must include the number of individuals of each species. The table with the plot attributes should include spatial coordinates of plots and information on experimental treatments and/or environmental variables. By default spatial coordinates are assumed to be named "x" and "y" for the easting and northing respectively unless otherwise specified by the argument `coord_names`. If the supplied coordinates are latitudes and longitudes be sure to set the argument `latlong = TRUE` which turns on great circle distances rather than Euclidean distances when computing the spatial, sample-based rarefaction (sSBR). If temporal trends are of interest rather than simply only supply one coordinate which represents time. # Invasion case study We will examine a case study on the effects of an invasive plant _Lonicera maackii_ on understory plant biodiversity in a Missouri woodland (Powell et al. 2003). This data was reanalyzed in McGlinn et al. (2018). This data set is included in `mobr`and available after loading the library ```{r load pkg and data, message=FALSE} library(mobr) library(dplyr) library(ggplot2) data(inv_comm) # Community matrix data(inv_plot_attr) # Plot attributes data.frame ``` ```{r examine data} str(inv_comm) head(inv_plot_attr) ``` The plot attributes include the information if a plot is located in invaded or uninvaded sites as well as the spatial xy coordinates of a plot. ## Preparing data In order to analyze the data with `mobr` the two data tables have to be combined into on single object ```{r make mob_in} inv_mob_in <- make_mob_in(inv_comm, inv_plot_attr, coord_names = c('x', 'y')) inv_mob_in ``` # Exploratory data analysis The package `mobr` offers functions for exploratory data analysis and visualization. First let's look at the spatial, sample-based rarefaction (sSBR) in which samples are added depending on their spatial proximity. ```{r, fig.width = 5, fig.height = 5} plot_rarefaction(inv_mob_in, 'group', ref_level = 'uninvaded', 'sSBR', lwd = 4) ``` We can see that invasion decreases richness and that the magnitude of the effect depends on scale. Let's dig in further to see if we can better understand exactly what components of the community are changing due to invasion. First we look at the individual rarefaction curves for each treatment which only reflect the shape of the SAD (i.e., no individual density or aggregation effects): ```{r, fig.width = 7, fig.height=4} oldpar <- par(no.readonly = TRUE) par(mfrow = c(1, 2)) plot_rarefaction(inv_mob_in, 'group', 'uninvaded', 'IBR', leg_loc = 'bottomright') par(oldpar) ``` Visually you can see there are some big differences in total numbers individuals (i.e., how far on the x-axis the curves go). We can also see that for small numbers of individuals the invaded sites are actually more diverse! This is a bit surprising and it implies that the invaded sites have greater evenness. We can directly examine the species abundance distribution (SAD): ```{r, fig.width = 7, fig.height=4} oldpar <- par(no.readonly = TRUE) par(mfrow = c(1, 2)) plot_abu(inv_mob_in, 'group', type = 'rad', scale = 'alpha', log = 'x') plot_abu(inv_mob_in, 'group', type = 'rad', scale = 'gamma' , log = 'x') par(oldpar) ``` The SADs suggest that the invaded site has greater evenness in its common species (i.e., flatter left hand-side of rank curve). # Two-scale analysis There are a myriad of biodiversity indices. We have attempted to chose a subset of metrics that can all be derived from the individual rarefaction curve which capture the aspects of biodiversity we are most interested in, namely: 1. Numbers of individuals (i.e., density effects) 2. The distribution of rarity and commonness (i.e., the SAD) 3. The spatial patchiness or aggregation of conspecifics. The metrics we have selected are: * **N - Number of individuals** * **S - Observed species richness** * **S_n - Rarefied species richness (Hurlbert 1971)** * **S_C - coverage based species richness (Chao & Jost 2012)** * **PIE - Probability of Interspecific Encounter (Hurlbert 1971)** * **S_PIE - Effective number of species based on PIE (Jost 2007)** Each of these metrics can be computed for either the sample or group scale individual rarefaction curves as shown in the figure below:
![](../vignettes/metrics.png){ width=100%}
The ratio of a given biodiversity metric at the group and sample scales (i.e., $\beta_S = S_{group} / S_{sample}$ can provide simple measures of species turnover or \(\beta\)-diversity. Depending on which metric the \(\beta\)-diversity is derived from will determine what aspects of species turnover it is most sensitive too (McGlinn et al. 2018, Chase et al. 2018). For the invasion dataset we will examine all of the metrics except for coverage based richness (*S_C*) because some of the samples in this dataset are empty which makes that calculation difficult for this example (see below for examples in which *S_C* is computed). ```{r compute diversity indices, eval = TRUE} indices <- c('N', 'S', 'S_n', 'S_PIE') inv_div <- tibble(inv_comm) %>% group_by(group = inv_plot_attr$group) %>% group_modify(~ calc_comm_div(.x, index = indices, effort = 5, extrapolate = TRUE)) ``` ```{r, echo = FALSE} #load('../vignettes/inv_stats.Rdata') ``` We can examine the `inv_div` object ```{r} head(inv_div) ``` There are also functions for plotting the indices at the sample and group levels. First let's examine species richness: ```{r, fig.width = 7, fig.height = 3.5} plot_comm_div(inv_div, 'S') ``` Invasion appears to decrease local sample diversity but not gamma diversity. Somewhat surprisingly it appears to increase \(\beta\)diversity. It is possible to compute beta diversity directly rather than using the `calc_comm_div` function by using the `calc_beta_div` function. Here is a quick demonstration. We can calculate the beta diversity for raw richness (`S`), rarefied richness at 5 individuals (`S_n`), and the richness transformed version of the probability of interspecific encounter (`S_PIE`). ```{r, beta demo} calc_beta_div(inv_comm, c('S', 'S_n', 'S_PIE'), effort = 5) ``` One of the major effects we observed in the individual rarefaction curve was that the invaded sites appeared to have fewer individuals, let's examine the test of that: ```{r, fig.width = 7, fig.height = 3.5} plot_comm_div(inv_div, 'N') ``` Clearly our intuition was correct there is a very strong negative effect of invasion on N. So it appears that the changes we observed in S may be due to the negative effect of invasion on N. Let's examine `S_n` to test this: ```{r, fig.width = 7, fig.height = 3.5} plot_comm_div(inv_div, 'S_n') ``` ```{r, fig.width = 7, fig.height = 3.5} plot_comm_div(inv_div, 'S_PIE') ``` We can also plot `f_0` but for this dataset this metric does not show strong patterns so we'll stop here for now. If you want to plot all the biodiversity metrics at once you can simply use: ```{r, eval=FALSE} plot_comm_div(inv_div) ``` # Multi-scale analysis The continuous scale analysis using `mobr` aims at disentangling the consequences of three biodiversity components on observed changes in species richness 1. **Species abundance distribution (SAD)** 2. **Number of individuals (N)** 3. **Aggregation (clumping) of conspecific individuals** To accomplish this we use three different rarefaction curves which each capture different aspects of community structure:
![](../vignettes/rarefaction_curves.png){ width=50%}
If we examine the difference between each of these curves in our two treatments we can learn how the treatment influences richness via its effects on different components of community structure.
![](../vignettes/deltaS_concept.png){ width=75%}
We can carry out this analysis in `mobr` using the function `get_delta_stats`. For the sake of speed we'll run the analysis with just 20 permutations but at least 200 are recommended for actual applications. ```{r multi-scale analysis, eval=FALSE} inv_mob_in <- make_mob_in(inv_comm, inv_plot_attr, coord_names = c('x', 'y')) inv_deltaS <- get_delta_stats(inv_mob_in, 'group', ref_level='uninvaded', type='discrete', log_scale=TRUE, n_perm = 199) ``` ```{r, echo=FALSE} load('../vignettes/inv_deltaS.Rdata') ``` The best way to examine the contents of this object is to plot it. First let's examine the three rarefaction curves: ```{r, fig.width=7, fig.height=3.5} plot(inv_deltaS, stat = 'b1', scale_by = 'indiv', display='S ~ effort') ``` Now let's consider the effect sizes as a function of scale ```{r, fig.width=7, fig.height=3.5} plot(inv_deltaS, stat = 'b1', scale_by = 'indiv', display='stat ~ effort') ``` The grey polygons above represent the 95% quantile for the null models of no treatment effect. # Herbaceous woodland plants response to fire Meyers et al. (2015) collected data on woody plants in the Missouri Ozarks. They collected 26 sites that had been recently burned and 26 sites that they considered unburned. We will reanalyze their data using the MoB framework. ```{r load fire data} # plant community in response to a prescribed fire treatment in a # central US woodland data(fire_comm) data(fire_plot_attr) fire_mob_in <- make_mob_in(fire_comm, fire_plot_attr, coord_names = c('x', 'y')) ``` Now that the data is loaded and we have created a `mob_in` object we can begin visualizing the patterns in the dataset. It is always important to graph any dataset before attempting to analyze its patterns. We will start our exploration with the sSBR. ```{r} par(mfrow=c(1,3)) plot_rarefaction(fire_mob_in, 'group', ref_level = 'unburned', 'IBR', leg_loc = NA) plot_rarefaction(fire_mob_in, 'group', ref_level = 'unburned', 'sSBR', leg_loc = 'bottomright') ``` The sSBR indicates that the unburned site has higher species richness across all scales and it appears that the two sites differ in the rate of diversity accumulation at certain scales. Let's examine the species-abundance distributions to see if they may have shifted and may be partially responsible for this observed shift in S. ```{r} oldpar <- par(no.readonly = TRUE) par(mfrow = c(1, 2)) plot_abu(fire_mob_in, 'group', ref_level = 'unburned', 'rad', leg_loc = 'topright',) plot_abu(fire_mob_in, 'group', ref_level = 'unburned', 'sad', leg_loc = NA) par(oldpar) ``` Above we plot the SAD in two different ways. The plot on the left is the standard rank curve and the plot on the right is the cumulative curve which McGill (2011) advocates for using. The cumulative curve makes it particularly clear that the shape (i.e., evenness) between the treatments is really similar, but in the rank curve it does appear that the burned site does some sites that are highly dominated by a single species. Let's use the two-scale analysis to continue this exploration of the data and to test some hypotheses regarding which components are driving this decrease in S in the burned sites. ```{r} indices <- c('N', 'S', 'S_C', 'S_n', 'S_PIE') fire_div <- tibble(fire_comm) %>% group_by(group = fire_plot_attr$group) %>% group_modify(~ calc_comm_div(.x, index = indices, effort = 5, extrapolate = TRUE)) ``` Let's graphically examine the metrics: ```{r, fig.width = 7, fig.height = 3.5} plot_comm_div(fire_div) ``` It is clear that the fire decreased abundance of woody plants - that is likely not a surprise. The question remains if this drove the decreases in *S* we observed in the sSBR and which are obvious in the second set of panels on *S*. It is worth noting that that the decrease in *S* is only significant at the \(\alpha\)-scale, and that fire actually increases \(\beta\)-diversity. When we examine the pattern of *Sn* we see that the effect size is half of what it was for *S* at the \(\alpha\)-scale. indicating that the decrease in *N* is helping to decrease *S*. Also the significant increase in \(\beta_S\) is not as strong in \(\beta_{Sn}\) which indicates: 1) that the burned sites have higher turnover in part because they have fewer individuals, and 2) that there is a detectable increase in intra-specific aggregation in the burned sites because \(\beta_{Sn}\) is still higher in the burned site. *S*PIE indicates that their is a shift towards higher evenness in the unburned sites that is stronger at the \(\gamma\)-scale. So from the two-scale analysis we can conclude that fire: 1) decreases *N* which decreases *S*, 2) increases aggregation, and 3) decreases evenness. Let's now look at the multi-scale analysis to examine the scale-dependence and relative strength of these effects. Note here we set the argument `log_scale` to `TRUE` because there is a large number of individuals and it is not necessary to compute the differences between the curves at all possible scales. ```{r, eval = FALSE} fire_deltaS <- get_delta_stats(fire_mob_in, 'group', ref_level = 'unburned', type = 'discrete', log_scale = TRUE, n_perm = 199, overall_p = TRUE) ``` ```{r, echo = FALSE} load('../vignettes/fire_deltaS.Rdata') ``` This time we'll examine the full suite of graphics for the multimetric analysis. ```{r, fig.width = 7, fig.height = 6, } plot(fire_deltaS, stat = 'b1', scale_by = 'indiv') ``` These confirm our expectations based on the two-scale analysis. The IBR indicates the increase in evenness in the unburned site particularly at coarse scales, the nsSBR illustrates that the N-effect is fairly strong, and the sSBR shows that their are some effects of aggregation. The delta effect curves can show this more clearly. These plots illustrate that: 1. None of these individual effects are of a really large magnitude ((\(\Delta\)*S* is never more or less than 2 species). 2. The *N*-effect is the largest driver of decreases in *S*. 3. The increased aggregation in the burned sites is only really relevant at fine spatial scales. 4. Although evenness is lower in the burned sites it is only outside of the acceptance intervals at the finest spatial scales. # Enrichment effect on aquatic experimental cattle tank communities Here we reanalyze data from Chase (2010) in which freshwater aquatic communities of macroinvertebrates and juvenile amphibians in artificial pools (i.e., cattle tanks) were enriched with phosphorus or not. ```{r} data(tank_comm) data(tank_plot_attr) tank_mob_in <- make_mob_in(tank_comm, tank_plot_attr, coord_names = c('x', 'y')) ``` The data is loaded let's examine the sSBR. ```{r, fig.width = 5, fig.height = 5} plot_rarefaction(tank_mob_in, 'group', ref_level = 'low', 'sSBR') ``` Clearly enrichment had a big effect on species richness as it is larger in the "high" productivity (i.e., enriched) community compared to the "low" community at all scales. It is rare to see such a strong asymptote on a rarefaction curve but it does appear that no new species were being encountered in the low community. ```{r, fig.width = 7, fig.height=4} oldpar <- par(no.readonly = TRUE) par(mfrow = c(1, 2)) plot_abu(tank_mob_in, 'group', ref_level = 'low', 'rad') plot_abu(tank_mob_in, 'group', ref_level = 'low', 'sad', leg_loc = NA) par(oldpar) ``` The SAD plots indicate only a modest shift in the evenness component of the SAD although clearly from the sSBR we know that the size of the species pool is larger in the high community. ```{r} indices <- c('N', 'S', 'S_C', 'S_n', 'S_PIE') tank_div <- tibble(tank_comm) %>% group_by(group = tank_plot_attr$group) %>% group_modify(~ calc_comm_div(.x, index = indices, effort = 5, extrapolate = TRUE)) ``` Let's graphically examine the metrics: ```{r, fig.width = 7, fig.height = 3.5} plot_comm_div(tank_div) ``` The two-scale analysis indicates that enrichment does not affect *N* but that it does increase *S* but only at the coarse scale. The higher productivity sites had higher \(\beta_S\) indicating that species turnover was higher due to enrichment. The analysis of *S**n* indicated that this increase in turnover was due in part to an increase in spatial aggregation in the enriched tanks. *S*PIE indicated that the pattern of evenness particularly in common species was not changed by enrichment (as we observed in the plots of the SADs). ```{r, eval = FALSE} tank_deltaS <- get_delta_stats(tank_mob_in, 'group', ref_level = 'low', inds = 10, log_scale = TRUE, type = 'discrete', n_perm=199) ``` ```{r, echo = FALSE} load('../vignettes/tank_deltaS.Rdata') ``` ```{r, fig.width = 7, fig.height = 6} plot(tank_deltaS, stat = 'b1') ``` The delta S plot indicates that enrichment increased intra-specific aggregation which resulted in decreased S, did not influence S via N, and did increase the size of the species pool which our framework identifies as an SAD effect. We know from our other analyses that this was not due to large changes in evenness. This example is particularly interesting because: 1) there is absolutely no density effect, 2) the aggregation effects and SAD effects cancel each other out so strongly that a standard analysis of plot *S* indicates no change due to enrichment when obviously there are large changes to the community influencing richness. # References Chase, J. M. (2010). Stochastic community assembly causes higher biodiversity in more productive environments. Science, 328: 1388-1391. Chao, A., and L. Jost. 2012. Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size. Ecology 93:2533–2547. Chiu, C.-H., Wang, Y.-T., Walther, B.A. & Chao, A. (2014). An improved nonparametric lower bound of species richness via a modified good-turing frequency formula. Biometrics, 70, 671-682. Gotelli, N.J. & Colwell, R.K. (2001). Quantifying biodiversity: procedures and pitfalls in the measurement and comparison of species richness. Ecology letters, 4, 379–391 Hurlbert, S.H. (1971). The Nonconcept of Species Diversity: A Critique and Alternative Parameters. Ecology, 52, 577–586 Jost, L. (2007). Partitioning Diversity into Independent Alpha and Beta Components. Ecology, 88, 2427-2439. Myers, J. A., Chase, J. M., Crandall, R. M., & Jiménez, I. (2015). Disturbance alters beta‐diversity but not the relative importance of community assembly mechanisms. Journal of Ecology, 103(5), 1291-1299. McGlinn, D.J. X. Xiao, F. May, S. Blowes, J. Chase, N. Gotelli, T. Knight, B. McGill, and O. Purschke. 2019. MoB (Measurement of Biodiversity): a method to separate the scale-dependent effects of species abundance distribution, density, and aggregation on diversity change. Methods in Ecology and Evolution. 10:258–269. https://doi.org/10.1111/2041-210X.13102 Powell, K.I., Chase, J.M. & Knight, T.M. (2013). Invasive Plants Have Scale-Dependent Effects on Diversity by Altering Species-Area Relationships. Science, 339, 316–318.