---
title: "Temporal Diversity Indices"
author: "Lauren M. Hallett, Dan F. B. Flynn, Sydney K. Jones"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
bibliography: biblio.bib
vignette: >
  %\VignetteIndexEntry{Temporal diversity indices}
  %\VignetteEngine{knitr::rmarkdown}
  %\SweaveUTF8
  %\VignetteEncoding{UTF-8}
  %\usepackage[utf8]{inputenc}
  
  
---

## 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 [@Collins2008; @cleland2013]. 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 [@Collins2008]. `codyn` includes three functions to characterize temporal shifts in species identity and rank abundances over time:

- `turnover` calculates total turnover as well as the proportion of species that either appear or disappear between time points.

- `rank_shift` quantifies 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_change` analyzes differences in species composition between samples at increasing time lags. It reflects the rate of directional change in community composition. The associated function `rate_change_interval` returns the full set of community distance values and associated time lag intervals 
to use in visualization.


```{r echo=FALSE}
library(knitr)
options(digits = 3)
```

```

## 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. 

```{r echo=FALSE}
library(codyn)
library(knitr)
data(collins08)
kable(head(collins08))
```

## 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 results='asis'}
KNZ_turnover <- turnover(df = collins08, 
                         time.var = "year", 
                         species.var = "species", 
                         abundance.var = "abundance", 
                         replicate.var = "replicate")
```
```{r echo=FALSE}
kable(head(KNZ_turnover))
```

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.

```{r results='asis', error = TRUE}
KNZ_turnover_agg <- turnover(df = collins08, 
                          species.var = "species",
                          time.var = "year",
                          abundance.var = "abundance")
```

### 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.
```{r results='asis'}
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. 
```{r results='asis'}
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.

```{r fig.width = 7, fig.height = 4, echo=FALSE}
library(ggplot2)
ggplot(KNZ_turnover, aes(x=year, y=total, color=replicate)) + geom_line(size = 2) + theme_bw()
```

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.

```{r fig.width = 7, fig.height = 4, echo=FALSE}
KNZ_appearance$metric<-"appearance"
names(KNZ_appearance)[1]="turnover"

KNZ_disappearance$metric<-"disappearance"
names(KNZ_disappearance)[1]="turnover"

KNZ_appdisapp<-rbind(KNZ_appearance, KNZ_disappearance)

ggplot(KNZ_appdisapp, aes(x=year, y=turnover, color=replicate)) + geom_line(size = 2) + theme_bw() + facet_wrap(~metric) + ggtitle("Species appearances and disappearances \n Annually burned vs unburned plots, Konza
          ")
```


## 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 [@Collins2008]. 

To illustrate, generate some sample data:

```{r}

## 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:
```{r, fig.width = 7, fig.height = 4}
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:
```{r, fig.width = 5, fig.height = 5}
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.

```{r, fig.width = 7, fig.height = 5}

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.

```{r}
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))
```
```{r echo=FALSE}
kable(head(KNZ_rankshift))
```

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:

```{r, error = TRUE}
KNZ_rankshift_agg <- rank_shift(df = collins08,
                        time.var = "year",
                        species.var = "species",
                        abundance.var = "abundance")
```

### 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. 
```{r fig.width = 7, fig.height = 4}
# 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 [@collins2000].


```{r}
    rate.res <- rate_change(collins08,  
                    time.var= "year", 
                    species.var= "species", 
                    abundance.var= "abundance", 
                    replicate.var = "replicate")
```

```{r, echo = F}
rownames(rate.res) = NULL
```

```{r, fig.width = 7, fig.height = 5}
kable(rate.res)

```

### 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.

```{r, fig.width = 7, fig.height = 5}

#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)

```

## Citations


