Temporal Diversity Indices

Lauren M. Hallett, Dan F. B. Flynn, Sydney K. Jones

2020-11-30

Overview

Many measures of community structure, such as diversity indices and rank-abundance curves, represent “snapshots in time” - which do not capture the temporal dynamics of ecological systems. For example, species richness (i.e., the number of species present) is a common metric to characterize ecological communities. However, species richness may be a poor indicator of community change if there is turnover in the identity but not the number of species present over time (Collins et al. 2008; Cleland et al. 2013). Similarly, species rank abundances can reorder through time, even when the pool of species present remains relatively stable. Within-year rank abundance curves fail to capture this dynamic, and other tools are needed to quantify within-community reordering without masking informative temporal shifts in internal dynamics (Collins et al. 2008). codyn includes three functions to characterize temporal shifts in species identity and rank abundances over time:


## Example dataset
`codyn` utilizes a subset of the data presented by Collins et al. [-@Collins2008]. The `collins08` dataset includes plant composition from the Konza Prairie Long-term Ecological Research site in Manhattan, KS. It spans 18 years and includes data collected from two locations, one that is annually burned and one that is unburned. 


|replicate       | year|species  | abundance|
|:---------------|----:|:--------|---------:|
|annually burned | 1984|achimill |     0.050|
|annually burned | 1984|ambrpsil |     1.350|
|annually burned | 1984|amorcane |     1.325|
|annually burned | 1984|andrgera |    76.000|
|annually burned | 1984|antenegl |     0.050|
|annually burned | 1984|apoccann |     0.075|

## Species turnover
Species turnover represents a temporal analog to species richness. It is described by the total turnover in species (appearances and disappearances) as well as the proportion of species that either appear or disappear between time points [@Collins2008; @cleland2013]. 

### Total turnover
The function `turnover` calculates three metrics of species turnover: total turnover, appearances, and disappearances.

The default metric `total` refers to total turnover, which calculates the proportion of species that differ between time points as:

$$ Total\; turnover = \frac{Species\; gained\; +\; Species\; lost}{Total\; species\; observed\; in\; both\; timepoints} $$

The metric was introduced by MacArthur and Wilson [-@MacArthur1963] and modified by Diamond [-@diamond1969] to be expressed as a proportion in order to compare turnover between sites that differ in species richness.

`turnover` requires a data frame with columns for species, time and abundance, and includes an optional argument to specify a column for spatial replicates.

```r
KNZ_turnover <- turnover(df = collins08, 
                         time.var = "year", 
                         species.var = "species", 
                         abundance.var = "abundance", 
                         replicate.var = "replicate")
total year replicate
0.237 1985 annually burned
0.302 1986 annually burned
0.242 1987 annually burned
0.230 1988 annually burned
0.185 1989 annually burned
0.224 1990 annually burned

Note that if the replicate.var column is not specified but the dataset in fact includes replicates, calling turnover will result in an error. For example, forgetting to specify the replicate.var, allowing it to assume the default value of NA, yields an error because multiple abundance values are included for the same species within a single year.

KNZ_turnover_agg <- turnover(df = collins08, 
                          species.var = "species",
                          time.var = "year",
                          abundance.var = "abundance")
## Error in check_single_onerep(df, time.var, species.var): Either data span multiple replicates with no replicate.var specified or multiple records within years for some species

Appearances and disappearances

Total turnover incorporates both species appearances and disappearances, but sometimes it is useful to parse their relative contribution. For example, a time point in which many species appear may reflect a different ecological story than a time point in which many species drop from the system, but the total turnover in both scenarios may be similar.

Specifying metric="appearance" will return the proportion of species that appeared relative to the total number of species observed in both time points. As before, spatial replicates can be specified with the replicate argument; setting replicate.var=NA will calculate across the full dataset.

KNZ_appearance <- turnover(df = collins08, 
                        time.var = "year",
                        species.var = "species",
                        abundance.var = "abundance",
                        replicate.var = "replicate",
                        metric = "appearance")

Similarly, specifying metric="disappearance" will return the proportion of species that disappeared relative to the total number of species observed in both time points.

KNZ_disappearance <- turnover(df = collins08,
                        time.var = "year",
                        species.var = "species",
                        abundance.var = "abundance",
                        replicate.var = "replicate",
                        metric = "disappearance")

Turnover at Konza

Total turnover indicates that in general there were greater fluctuations in the species present in the unburned than burned location at Konza.

Parsing total turnover by appearances and disappearances indicates that there were generally a higher proportion of appearances in the unburned than burned location. In addition, the strong divergence in total turnover between locations in 1995 was because the unburned location lost a high proportion of species that year, while the annually burned locations lost none.

Species rank shifts

Rank clocks: visualizing species rank shifts

Rank clocks are a useful way to visualize the degree to which species reorder over time. They plot the rank order of abundance of each species over time in a circle, starting with a vertical axis at 12 o’clock (Collins et al. 2008).

To illustrate, generate some sample data:

## Generate some sample data 
yr = 1977:2003
sp1 = .3*sin(yr) + rnorm(length(yr), 0, .1) + -.05*yr + 150
sp2 = .2*sin(yr) + rnorm(length(yr), 0, .1) + -.01*yr + 70
sp3 = .2*sin(yr) + rnorm(length(yr), 0, .1) + .01*yr + 30

dat = data.frame(year = rep(yr, 3), 
           species = gl(3, length(yr), labels = c("sp1","sp2","sp3")),
           abundance = c(sp1, sp2, sp3))

A traditional, linear depiction of species rank shifts can be cluttered and visually confusing:

ggplot(dat, aes(year, abundance, color = species)) + 
  geom_line(size = 2) + theme_bw()

In contrast, a rank clock visualization highlights differences in the stability of species:

ggplot(dat, aes(year, abundance, color = species)) + 
  geom_line(size = 2) + coord_polar() + theme_bw()

Rank clocks at Konza

Rank clocks highlight that there has been tremendous reordering in the relative abundance of dominant species in the annually burned but not the unburned location at Konza. For example, big bluestem (Andropogon gerardii) decreased substantially in the annually burned plot over time but remained stable and dominant in the unburned plot.

aggdat <- aggregate(abundance ~ species * year * replicate, 
                    data = subset(collins08, 
                                    species == "andrgera" |
                                    species == "andrscop" | 
                                    species == "poaprat"| 
                                    species == "sorgnuta"), 
                    FUN = mean)

ggplot(aggdat, aes(year, abundance, color = species)) + 
  geom_line(size = 2) + coord_polar() + theme_bw() + facet_wrap(~replicate) +
  ggtitle("Dominant species abundances \n Annually burned vs unburned plots, Konza \n")

Mean rank shifts

The rank_shift function describes relative changes in species rank abundances, which indicate the degree of species reordering between two time points. This metric is calculated as:

\[ MRS = {\sum_{i=1}^{N} (|R_i,t+1 - R_i,t|})/N \]

where \(N\) is the number of species in common in both time points, \(t\) is the time point, \(R_i,t\) is the relative rank of species \(i\) in time \(t\).

rank_shift requires a data frame with columns for species, time and abundance, and includes an optional argument to specify a column for spatial replicates.

KNZ_rankshift <- rank_shift(df=collins08,
                        time.var = "year",
                        species.var = "species",
                        abundance.var = "abundance", 
                        replicate.var = "replicate")
#Select the final time point from the returned time.var_pair
KNZ_rankshift$year <- as.numeric(substr(KNZ_rankshift$year_pair, 6,9))
year_pair MRS replicate year
1984-1985 3.53 annually burned 1985
1985-1986 4.73 annually burned 1986
1986-1987 4.62 annually burned 1987
1987-1988 5.32 annually burned 1988
1988-1989 5.27 annually burned 1989
1989-1990 4.67 annually burned 1990

Calling rank_shift without specifying a replicate.var for a dataset that does in fact include replicates will result in an error. In the collins08 dataset, for example, if replicate.var = "replicate" is not given:

KNZ_rankshift_agg <- rank_shift(df = collins08,
                        time.var = "year",
                        species.var = "species",
                        abundance.var = "abundance")
## Error in check_single_onerep(df, time.var, species.var): Either data span multiple replicates with no replicate.var specified or multiple records within years for some species

Rank shifts at Konza

Calculating mean rank shifts highlights that the stability of communities diverged at Konza around 1992, with fewer rank shifts between species in the burned relative to the unburned area.

# Create a column with the final year from the returned time.var_pair
KNZ_rankshift$year <- as.numeric(substr(KNZ_rankshift$year_pair, 6, 9))

# Plot it
ggplot(KNZ_rankshift, aes(year, MRS, color=replicate)) + 
  geom_line(size= 2) + theme_bw()

Rate of community change

The rate_change function assesses the rate and pattern of variability within a community, which indicates whether species reordering over time is resulting in directional change. This function calculates differences in species composition between samples at increasing time intervals. Differences in species composition are characterized by Euclidean distances, which are calculated on pair-wise communities across the entire time series. For example, a data set with 6 time intervals will have distance values for five one-interval time lags (e.g., time 1 vs time 2, time 2 vs time 3 …), 4 two-interval time lags (e.g., time 1 vs time 3, time 2 vs time 4 …) and so forth. These distance values are regressed against the time lag interval. The slope of the regression line is an indication of the rate and direction of compositional change in the community (Collins, Micheli, and Hartt 2000).

    rate.res <- rate_change(collins08,  
                    time.var= "year", 
                    species.var= "species", 
                    abundance.var= "abundance", 
                    replicate.var = "replicate")
kable(rate.res)
replicate rate_change
annually burned 2.810
unburned 0.147

Rate change at Konza

The annually burned grassland location shown here has relatively greater directional change through time than the unburned location. This reflects the pattern highlighted in the rank clocks, in which, for example, the dominant species Andropogon gerardii decreased substantially over time in the annually burned location but remained stable in the unburned location.

#Use the rate_change_interval function to generate the full data frame of distances by time lag intervals
comm.res <- rate_change_interval(collins08, 
                              time.var = "year",
                              species.var = "species",
                              abundance.var = "abundance",
                              replicate.var = "replicate")

ggplot(comm.res, aes(interval, distance, color = replicate)) + facet_wrap(~replicate) + 
  geom_point() + theme_bw() + stat_smooth(method = "lm", se = F, size = 2)
## `geom_smooth()` using formula 'y ~ x'

Citations

Cleland, Elsa E., Scott L. Collins, Timothy L. Dickson, Emily C. Farrer, Katherine L. Gross, Laureano A. Gherardi, Lauren M. Hallett, et al. 2013. “Sensitivity of Grassland Plant Community Composition to Spatial Vs. Temporal Variation in Precipitation.” Ecology 94 (8): 1687–96. https://doi.org/10.1890/12-1006.1.

Collins, Scott L., Fiorenza Micheli, and Laura Hartt. 2000. “A Method to Determine Rates and Patterns of Variability in Ecological Communities.” Oikos 91 (2): 285–93. https://doi.org/10.1034/j.1600-0706.2000.910209.x.

Collins, Scott L., Katharine N. Suding, Elsa E. Cleland, Michael Batty, Steven C. Pennings, Katherine L. Gross, James B. Grace, Laura Gough, Joe E. Fargione, and Christopher M. Clark. 2008. “Rank Clocks and Plant Community Dynamics.” Article. Ecology 89 (12): 3534–41. https://doi.org/10.1890/07-1646.1.