Temporal Diversity Indices
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:
turnovercalculates total turnover as well as the proportion of species that either appear or disappear between time points.rank_shiftquantifies relative changes in species rank abundances by taking the sum difference of species ranks in consecutive time points. This metric goes hand-in-hand with “rank clocks,” a useful visualization tool of shifts in species ranks.rate_changeanalyzes differences in species composition between samples at increasing time lags. It reflects the rate of directional change in community composition. The associated functionrate_change_intervalreturns the full set of community distance values and associated time lag intervals to use in visualization.
## 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 speciesAppearances 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.
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.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
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:
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, Ri, 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 speciesRank 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")| 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'