diversity.Rmd

---
title: "diversity"
author: "Antonio J Perez-Luque"
date: "2021-05-03"
output: workflowr::wflow_html
editor_options:
  chunk_output_type: console
---

# Diversity analysis 

- Compute Simpson and Shannon diversity index. See diversity functions in vegan pkg. 

```{r, set-global-options}
knitr::opts_chunk$set(echo = TRUE, 
                      warning = FALSE, 
                      message = FALSE)
```

```{r load-pkgs}
library("tidyverse")
library("here")
library("vegan")
```

```{r}
passerine <- read_csv(here::here("data/passerine.csv")) 
passerine_abbreviations <- read_csv(here::here("data/passerine_abbrev.csv")) 


passerine.ab <- passerine %>% 
  rename(sp = sp.abb) %>% 
  dplyr::select(-especie) %>% 
  group_by(sp, year, habitat) %>% 
  summarise(ab_avg = round(mean(den, na.rm=TRUE),2), 
            sd = sd(den, na.rm = TRUE), 
            se = sd/sqrt(length(den)), 
            n = length(den))

df <- passerine.ab %>% 
  unite("site", habitat:year, sep="_", remove = FALSE) %>% 
  dplyr::select(-sd, -se, -n) %>% 
  pivot_wider(names_from = sp, 
              values_from = ab_avg, 
              values_fill = 0)
```


Create a function to compute Simpson and Shannon diversity indices

```{r, computeDiversity}
computeDiversity <- function(d, hab, ...){
  
  df.habitat <- d %>% 
    filter(habitat == hab) %>% 
  dplyr::select(-habitat, -site) %>% 
  dplyr::select_if(
    negate(function(x) is.numeric(x) && sum(x) == 0)) %>%
    column_to_rownames("year")
  

  shannon <- diversity(df.habitat, "shannon") %>%
    as.data.frame() %>% rownames_to_column("year") %>%
    rename(shannon = starts_with("."))

  simpson <- diversity(df.habitat, "simpson") %>%
    as.data.frame() %>%
    rownames_to_column("year") %>%
    rename(simpson = starts_with("."))

  diver <- shannon %>%
  inner_join(simpson) %>%
  mutate(habitat = hab)
  
  return(diver)
}
```


- Compute diversity index for each site 

```{r}
diver.summits <- computeDiversity(df, "cumbres")
diver.juniper <- computeDiversity(df, "enebral")
diver.robledal <- computeDiversity(df, "robledal")

diversity_index <- bind_rows(diver.summits, diver.juniper, diver.robledal)
write_csv(diversity_index, here::here("output/diversity/alfa_diversity.csv"))
```

```{r, plot-diversisty}
diversity_plot <- diversity_index %>% 
  filter(!(habitat == "cumbres" & year == 2013)) %>% 
  ggplot(
  aes(x=year, y=shannon, colour=habitat, group=habitat)) +
    geom_point() + geom_line() +
  geom_smooth(method = "lm", 
              se=FALSE, size = .5, linetype="dashed") + 
  ylab("Shannon Index (H')") + 
  theme_bw() + 
  theme(
    panel.grid = element_blank(),
    #panel.grid.major.x = element_blank(),
    strip.background = element_blank()
  )


ggsave(here::here("output/diversity/diversity_habitat.pdf"),
       width = 18, height = 12, units = "cm")
diversity_plot
dev.off()
```

# Compare three temporal times

- Generate three temporal times: 1983; 2009: mean for 2008 - 2009; and 2018: mean for 2017 - 2018. 

```{r, fig.width=3}
shannon <- diversity_index %>% 
  filter(year %in% c(1981:1985, 2008:2009, 2017:2018)) %>% 
  mutate(yearF = case_when(
    year %in% c(1981:1985) ~ as.character(1983), 
    year %in% c(2008:2009) ~ as.character(2009),
    year %in% c(2017:2018) ~ as.character(2018)
  )) %>% 
  group_by(yearF, habitat) %>% 
  summarise(mean = mean(shannon, na.rm=TRUE), 
            sd = sd(shannon, na.rm=TRUE),
            se = sd/sqrt(length(shannon)))

shannon %>% 
  ggplot(aes(x=yearF, y=mean, group=1)) + 
  geom_point() +
  geom_line() + 
  # geom_bar(stat="identity", fill="black") + 
  geom_errorbar(aes(ymax = mean + se, ymin= mean - se), 
                width = .2, size = .3) + 
  facet_wrap(~habitat, ncol = 1) + 
  theme_bw() +
  theme(
    panel.grid.minor = element_blank(),
    strip.background = element_blank()
  ) + 
  ylab("Shannon Index (H')") + xlab("")
```

## Richness 

```{r}
richness <- passerine.ab %>% 
  group_by(habitat, year) %>% 
  summarise(richness = length(sp))
write_csv(richness, here::here("output/diversity/richness.csv"))
```

```{r}
richness %>% 
  ggplot(aes(as.factor(year), richness)) +
  geom_point() + 
  geom_path(aes(group=habitat)) + 
  facet_wrap(~habitat, nrow = 3) + 
  ylab("Richness (n sps)") + xlab("") + 
  theme_bw() + 
  theme(
    panel.grid.minor = element_blank(),
    strip.background = element_blank()
  )
```


```{r}
richness %>% 
  filter(year %in% c(1981:1985, 2008:2009, 2017:2018)) %>% 
  mutate(yearF = case_when(
    year %in% c(1981:1985) ~ as.character(1983), 
    year %in% c(2008:2009) ~ as.character(2009),
    year %in% c(2017:2018) ~ as.character(2018)
  )) %>% 
  group_by(yearF, habitat) %>% 
  summarise(mean = mean(richness, na.rm=TRUE), 
            sd = sd(richness, na.rm=TRUE),
            se = sd/sqrt(length(richness))) %>% 
  ggplot(aes(x=yearF, y=mean, group=1)) + 
  geom_point() +
  geom_line() + 
  # geom_bar(stat="identity", fill="black") + 
  geom_errorbar(aes(ymax = mean + se, ymin= mean - se), 
                width = .2, size = .3) + 
  facet_wrap(~habitat, ncol = 1) + 
  theme_bw() +
  theme(
    panel.grid.minor = element_blank(),
    strip.background = element_blank()
  ) 

```