## ----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 |>
  select(TAXON, GENERIC, SPECIFIC)

## -------------------------------------------------------------------------------------------------
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", seq_len(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 |>
  pivot_longer(cols = starts_with("POOLRUN"), names_to = "POOLRUN",
               values_to = "POOLID") |>
  mutate(POOLRUN = as.integer(parse_number(POOLRUN))) |>
  filter(!is.na(POOLID)) |>
  select(POOLRUN, POOLID, HABITAT, AREA, YEAR, ID, TAXON, COUNT)
oosterschelde

## -------------------------------------------------------------------------------------------------
d <- oosterschelde |>
  group_by(HABITAT, YEAR, POOLRUN, POOLID) |>
  distinct(ID, AREA) |>
  summarise(AREA = sum(AREA), .groups = "drop")
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", linewidth = 1,
             alpha = 0.5) +
  geom_linerange(
    mapping = aes(x = AREA, ymin = 0, ymax = n),
    colour = "seashell4", linewidth = 2
  ) +
  scale_x_continuous(name = expression("area pooled sample (m"^2 * ")"))

## -------------------------------------------------------------------------------------------------
d <- oosterschelde |>
  group_by(HABITAT, YEAR, POOLRUN, POOLID) |>
  summarise(S = species_richness(taxon = TAXON, count = COUNT),
            .groups = "drop")
d

## -------------------------------------------------------------------------------------------------
d <- d |>
  group_by(HABITAT, YEAR) |>
  summarise(S = mean(S), .groups = "drop")
d

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

## -------------------------------------------------------------------------------------------------
total_abundance(count = d$COUNT)

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

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

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

## -------------------------------------------------------------------------------------------------
rygg(taxon = d$TAXON, count = d$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)
  ) |>
  pivot_longer(cols = all_of(c("margalef", "SN_rygg", "SN_adj")),
               names_to = "index")
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_transform(x = "log10") +
  facet_grid(index ~ ., scales = "free_y")

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

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

## ----echo=FALSE, fig.width=7, fig.height=4, out.width=500, echo=FALSE-----------------------------
n <- seq_len(total_abundance(count = d$COUNT))
ESn <- sapply(X = n,
              FUN = \(n) hurlbert(taxon = d$TAXON, count = d$COUNT, n = n))
ggplot(data = tibble(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])))

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

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

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

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

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

## ----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 = \(a) hill(taxon = d$TAXON, count = d$COUNT, a = a))
ggplot(data = tibble(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))

## -------------------------------------------------------------------------------------------------
ambi(taxon = d$TAXON, count = d$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()

## -------------------------------------------------------------------------------------------------
iti(taxon = d$TAXON, count = d$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),
    .groups = "drop"
  )

## -------------------------------------------------------------------------------------------------
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) |>
  pivot_longer(cols = starts_with("POOLRUN"), names_to = "POOLRUN",
               values_to = "POOLID") |>
  mutate(POOLRUN = as.integer(parse_number(POOLRUN))) |>
  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),
            .groups = "drop") |>
  group_by(HABITAT, YEAR, POOLRUN) |>
  summarise(S = mean(S), .groups = "drop_last") |>
  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)