benthos.R
In benthos: Marine Benthic Ecosystem Analysis

## ----setup, include = FALSE-----------------------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = NA_character_
)
options(dplyr.width = 100, width = 100)
# see e-mail Kurt Hornik 2019-03-05T12:59
suppressWarnings(RNGversion("3.5.0"))
set.seed(314)

## -------------------------------------------------------------------------------------------------
library(benthos)

## ---- message=FALSE-------------------------------------------------------------------------------
library(dplyr)
library(tidyr)
library(readr)
library(ggplot2)

## -------------------------------------------------------------------------------------------------
data(oosterschelde)

## -------------------------------------------------------------------------------------------------
oosterschelde

## ---- eval=FALSE----------------------------------------------------------------------------------
#  ?oosterschelde

## -------------------------------------------------------------------------------------------------
oosterschelde <- oosterschelde %>%
    filter(months(DATE) %in% c("August", "September"))

## -------------------------------------------------------------------------------------------------
oosterschelde %>% 
    mutate(COMPLIANT = is_accepted(taxon = TAXON)) %>%
    select(OBJECTID, HABITAT, DATE, TAXON, COUNT, COMPLIANT)

## -------------------------------------------------------------------------------------------------
oosterschelde %>% 
    mutate(COMPLIANT = is_accepted(taxon = TAXON)) %>%
    summarise(
        N_RECORDS = n(),
        N_COMPLIANT = sum(COMPLIANT),
        N_MISSING = N_RECORDS - N_COMPLIANT
    )

## -------------------------------------------------------------------------------------------------
oosterschelde %>% 
    filter(!is_accepted(taxon = TAXON))
oosterschelde <- oosterschelde %>% 
    filter(is_accepted(taxon = TAXON))

## -------------------------------------------------------------------------------------------------
is_accepted(taxon = "Petricola pholadiformis")
as_accepted(taxon = "Petricola pholadiformis")
is_accepted(taxon = "Petricolaria pholadiformis")

## -------------------------------------------------------------------------------------------------
oosterschelde <- oosterschelde %>%
    mutate(TAXON = as_accepted(TAXON))

## -------------------------------------------------------------------------------------------------
oosterschelde %>% 
    filter(!is_accepted(taxon = TAXON)) %>%
    nrow

## -------------------------------------------------------------------------------------------------
oosterschelde <- oosterschelde %>%
    mutate(
        GENERIC  = generic_name(TAXON),
        SPECIFIC = specific_name(TAXON)
    )
oosterschelde

## -------------------------------------------------------------------------------------------------
is_binomen("Nephtys")
is_binomen("Nephtys cirrosa")

## -------------------------------------------------------------------------------------------------
oosterschelde <- oosterschelde %>%
    mutate(
        IS_GENUS = is.na(GENERIC),
        GENERIC = ifelse(IS_GENUS, TAXON, GENERIC)
    )

## -------------------------------------------------------------------------------------------------
oosterschelde %>%
    filter(IS_GENUS) %>%
    nrow

## -------------------------------------------------------------------------------------------------
oosterschelde <- oosterschelde %>%
    group_by(ID, GENERIC) %>%
    mutate(NEWCOUNT = genus_to_species(is_genus = IS_GENUS, count = COUNT)) %>%
    ungroup

## -------------------------------------------------------------------------------------------------
oosterschelde %>%
    filter(GENERIC == "Corophium", ID == 130) %>%
    arrange(TAXON)

## -------------------------------------------------------------------------------------------------
oosterschelde %>%
    summarise(sum(COUNT), sum(NEWCOUNT))

## -------------------------------------------------------------------------------------------------
oosterschelde <- oosterschelde %>%
    mutate(COUNT = NEWCOUNT) %>%
    select(ID, HABITAT, AREA, DATE, TAXON, COUNT) %>%
    filter(COUNT > 0)
oosterschelde

## -------------------------------------------------------------------------------------------------
d <- oosterschelde %>%
    group_by(AREA) %>%
    summarise(n = n())
d

## -------------------------------------------------------------------------------------------------
oosterschelde <- oosterschelde %>%
    mutate(YEAR = DATE %>% format("%Y"))

## -------------------------------------------------------------------------------------------------
n_pool_runs <- 10
pool_runs <- replicate(
    n = n_pool_runs, {
    oosterschelde %>% 
        group_by(HABITAT, YEAR) %>%
        mutate(
            POOLID = pool(
                sample_id = ID, 
                area = AREA, 
                target_area = c(0.09, 0.11)
            )
        ) %>%
        ungroup %>%
        select(POOLID)
    }
)

## -------------------------------------------------------------------------------------------------
names(pool_runs) <- paste0("POOLRUN", 1:n_pool_runs)
pool_runs <- pool_runs %>% as_tibble
pool_runs

## -------------------------------------------------------------------------------------------------
oosterschelde_orig <- oosterschelde
oosterschelde <- oosterschelde %>%
    bind_cols(pool_runs) %>% 
    as_tibble 
oosterschelde

## -------------------------------------------------------------------------------------------------
oosterschelde <- oosterschelde %>% 
    gather(key = "POOLRUN", value = "POOLID", starts_with("POOLRUN")) %>%
    mutate(POOLRUN = parse_number(POOLRUN) %>% as.integer) %>%
    filter(!is.na(POOLID)) %>%
    select(POOLRUN, POOLID, HABITAT, AREA, YEAR, ID, TAXON, COUNT)
oosterschelde

## -------------------------------------------------------------------------------------------------
d <- oosterschelde %>%
    group_by(HABITAT, YEAR, POOLRUN, POOLID) %>%
    select(ID, AREA) %>%
    distinct(ID, AREA) %>%
    summarise(AREA = sum(AREA))
d

## -------------------------------------------------------------------------------------------------
d <- d %>%
    select(AREA) %>%
    group_by(AREA) %>%
    summarise(n = n()) %>%
    arrange(AREA)
d

## ---- fig.retina=NULL, fig.width=5, fig.height=3, out.width=500, echo=FALSE-----------------------
ggplot(data = d) +
    geom_vline(xintercept = c(0.09, 0.11), colour = "red", size = 1, alpha = 0.5) +
    geom_linerange(
        mapping = aes(x = AREA, ymin= 0, ymax = n),
        colour = "seashell4", size = 2
    ) +
    scale_x_continuous(name = "area pooled sample (m2)")

## -------------------------------------------------------------------------------------------------
d <- oosterschelde %>%
    group_by(HABITAT, YEAR, POOLRUN, POOLID) %>%
    summarise(S = species_richness(taxon = TAXON, count = COUNT))
d

## -------------------------------------------------------------------------------------------------
d <- d %>%
    group_by(HABITAT, YEAR) %>%
    summarise(S = mean(S))
d

## -------------------------------------------------------------------------------------------------
d <- oosterschelde %>% 
    filter(HABITAT == "Polyhaline-Subtidal", YEAR == 2010, POOLRUN == 1, POOLID == 1) %>%
    select(TAXON, COUNT) %>%
    arrange(TAXON)
d

## -------------------------------------------------------------------------------------------------
d %>% total_abundance(count = COUNT)

## -------------------------------------------------------------------------------------------------
d %>% abundance(taxon = TAXON, count = COUNT) %>% as.matrix

## -------------------------------------------------------------------------------------------------
d %>% species_richness(taxon = TAXON, count = COUNT)

## -------------------------------------------------------------------------------------------------
d %>% margalef(taxon = TAXON, count = COUNT)

## -------------------------------------------------------------------------------------------------
d %>% rygg(taxon = TAXON, count = COUNT)

## ---- echo=FALSE, message=FALSE-------------------------------------------------------------------
x <- bind_rows(
    tibble(S =  1, N = c(1:9, 1:10*10)),
    tibble(S =  2, N = c(2:9, 1:10*10)),
    tibble(S =  3, N = c(3:9, 1:10*10)),
    tibble(S =  4, N = c(4:9, 1:10*10)),
    tibble(S =  5, N = c(5:9, 1:10*10)),
    tibble(S = 10, N = 1:10*10),
    tibble(S = 25, N = c(25, 3:10*10))) %>% 
    mutate(
        margalef = (S-1) / log(N),
        SN_rygg = log(S) / log(log(N)),
        SN_adj  = log(S) / log1p(log1p(N)),
        S = factor(S, ordered = TRUE)
    ) %>% 
    gather(key = "index", value = "value", margalef, SN_rygg, SN_adj)
x$value[!is.finite(x$value) | is.nan(x$value)] <- NA_real_

## ---- fig.width=7, fig.height=7, out.width=500, echo=FALSE, warning=FALSE-------------------------
ggplot(
    data = x, 
    mapping = aes(x = N, y = value, group = S, colour = S)
) +
    geom_path(size = 1, na.rm = TRUE) +
    geom_point(size = 3, na.rm = TRUE) +
    scale_x_continuous(breaks = c(1, 2, 3, 5, 10, 25, 50, 75, 100)) +
    coord_trans(x = "log10") +
    facet_grid(index~., scales = "free_y")

## -------------------------------------------------------------------------------------------------
d %>% rygg(taxon = TAXON, count = COUNT, adjusted = TRUE)

## -------------------------------------------------------------------------------------------------
d %>% hurlbert(taxon = TAXON, count = COUNT, n = 100)

## ---- echo=FALSE, fig.width=7, fig.height=4, out.width=500, echo=FALSE----------------------------
n <- seq_len(d %>% total_abundance(count = COUNT))
ESn <- sapply(X = n, FUN = function(n) {
    d %>% hurlbert(taxon = TAXON, count = COUNT, n = n)
})
ggplot(data = data.frame(n = n, ESn = ESn)) +
    geom_point(mapping = aes(x = n, y = ESn)) +
    scale_x_continuous(name = expression(italic(n))) +
    scale_y_continuous(name = expression(E(italic(S)[n])))

## -------------------------------------------------------------------------------------------------
d %>% simpson(taxon = TAXON, count = COUNT)

## -------------------------------------------------------------------------------------------------
d %>% hpie(taxon = TAXON, count = COUNT)

## -------------------------------------------------------------------------------------------------
1 - d %>% simpson(taxon = TAXON, count = COUNT)

## -------------------------------------------------------------------------------------------------
d %>% shannon(taxon = TAXON, count = COUNT)

## -------------------------------------------------------------------------------------------------
d %>% hill(taxon = TAXON, count = COUNT, a = 0)
d %>% hill(taxon = TAXON, count = COUNT, a = 1)
d %>% hill(taxon = TAXON, count = COUNT, a = 2)

## -------------------------------------------------------------------------------------------------
d %>% hill0(taxon = TAXON, count = COUNT)
d %>% hill1(taxon = TAXON, count = COUNT)
d %>% hill2(taxon = TAXON, count = COUNT)

## ---- fig.retina=NULL, fig.width=6, fig.height=4, out.width=500, echo=FALSE, warning=FALSE, message=FALSE----
a <- seq(from = 0, to = 2, by = 0.1)
N_a <- sapply(X = a, FUN = function(a) {
    d %>% hill(taxon = TAXON, count = COUNT, a = a)
})
ggplot(data = data.frame(a, N_a)) +
    geom_path(mapping = aes(x = a, y = N_a)) +
    geom_text(x = 0, y = 1, label = "<- rare species",   hjust = 0, vjust = 1, colour = "blue") +
    geom_text(x = 2, y = 1, label = "common species ->", hjust = 1, vjust = 1, colour = "blue") +
    scale_x_continuous(name = expression(italic(a))) +
    scale_y_continuous(name = expression(italic(N[a])), limits = c(0, NA))

## -------------------------------------------------------------------------------------------------
d %>% 
    ambi(taxon = TAXON, count = COUNT)

## -------------------------------------------------------------------------------------------------
d %>%
    mutate(HAS_GROUP = has_ambi(taxon = TAXON))

## -------------------------------------------------------------------------------------------------
d %>%
    mutate(HAS_GROUP = has_ambi(taxon = TAXON)) %>%
    summarise(percentage = 100 * sum(COUNT[!HAS_GROUP]) / sum(COUNT)) %>%
    as.numeric

## -------------------------------------------------------------------------------------------------
d %>% 
    iti(taxon = TAXON, count = COUNT)

## -------------------------------------------------------------------------------------------------
d %>%
    mutate(HAS_GROUP = has_iti(taxon = TAXON))

## -------------------------------------------------------------------------------------------------
d %>%
    mutate(HAS_GROUP = has_iti(taxon = TAXON)) %>%
    summarise(percentage = 100 * sum(COUNT[!HAS_GROUP]) / sum(COUNT)) %>%
    as.numeric

## -------------------------------------------------------------------------------------------------
oosterschelde %>% 
    group_by(HABITAT, YEAR, POOLRUN, POOLID) %>% 
    summarise(
        N = total_abundance(count = COUNT),
        S = species_richness(taxon = TAXON, count = COUNT),
        D = margalef(taxon = TAXON, count = COUNT),
        SN = rygg(taxon = TAXON, count = COUNT),
        SNa = rygg(taxon = TAXON, count = COUNT, adjusted = TRUE),
        H = shannon(taxon = TAXON, count = COUNT)
    )

## -------------------------------------------------------------------------------------------------
oosterschelde %>% 
    group_by(HABITAT, YEAR, POOLRUN, POOLID) %>% 
    summarise(
        N = total_abundance(., COUNT),
        S = species_richness(., TAXON, COUNT),
        D = margalef(., TAXON, COUNT),
        SN = rygg(., TAXON, COUNT),
        SNa = rygg(., TAXON, COUNT, adjusted = TRUE),
        H = shannon(., TAXON, COUNT)
    )

## -------------------------------------------------------------------------------------------------
oosterschelde <- oosterschelde_orig
n_pool_runs <- 100
pool_runs <- replicate(
    n = n_pool_runs, {
    oosterschelde %>% 
        group_by(HABITAT, YEAR) %>%
        mutate(
            POOLID = pool(
                sample_id = ID, 
                area = AREA, 
                target_area = c(0.09, 0.11)
            )
        ) %>%
        ungroup %>%
        select(POOLID)
    }
)
names(pool_runs) <- paste0("POOLRUN", 1:n_pool_runs)

d <- pool_runs %>% 
    as_tibble %>%
    bind_cols(oosterschelde) %>% 
    as_tibble %>%
    gather(key = "POOLRUN", value = "POOLID", starts_with("POOLRUN")) %>%
    mutate(POOLRUN = parse_number(POOLRUN) %>% as.integer) %>%
    filter(!is.na(POOLID)) %>%
    select(POOLRUN, POOLID, HABITAT, AREA, YEAR, ID, TAXON, COUNT) %>%
    group_by(HABITAT, YEAR, POOLRUN, POOLID) %>%
    summarise(S = species_richness(taxon = TAXON, count = COUNT)) %>%
    group_by(HABITAT, YEAR, POOLRUN) %>%
    summarise(S = mean(S)) %>%
    mutate(S_rm = cummean(S)) 
d

## ---- fig.retina=NULL, fig.width=7, fig.height=4, out.width=650, echo=FALSE, warning=FALSE--------
ggplot(data = d) +
    geom_point(mapping = aes(x = POOLRUN, y = S), col = "blue") +
    geom_path(mapping = aes(x = POOLRUN, y = S_rm), col = "red") +
    facet_grid(HABITAT ~ YEAR)