## ----echo=FALSE---------------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
options(rmarkdown.html_vignette.check_title = FALSE)
set.seed(123)

## ----setup, include=FALSE-----------------------------------------------------
library(estar)
library(tidyr)
library(dplyr)
library(ggplot2)
library(cowplot)
library(viridis)
library(forcats)
library(kableExtra)
library(vegan)
source("custom_aesthetics.R")

label_intensity <- function(intensity, mu) {
  paste0("Conc = ", intensity, " micro g/L")
}

## ----echo = FALSE, message = FALSE,  fig.height = 6, fig.width = 8, fig.cap = "Organism mean count (log axis) of the five functional groups over the course of ecotoxicological experiment. The vertical dashed line identifies the time when the insecticide was applied (at week 0), and facets represent the different concentrations (micro g/L)."----
(
  aquacomm_resps.logplot <- aquacomm_fgps |>
    tidyr::pivot_longer(
      cols = c(herb, detr_herb, carn, omni, detr),
      names_to = "gp",
      values_to = "abund"
    ) |>
    ## summarize to plot
    dplyr::group_by(time, treat, gp) |>
    dplyr::summarize(abund_logmean = log(mean(abund))) |>
    dplyr::ungroup() |>
    dplyr::mutate(
      gp = forcats::fct_recode(
        factor(gp),
        Herbivore = "herb",
        Detr_Herbivore = "detr_herb",
        Carnivore = "carn",
        Omnivore = "omni",
        Detritivore = "detr"
      )
    ) |>
    ggplot(aes(
      x = time,
      y = abund_logmean,
      group = gp,
      colour = gp
    )) +
    geom_line() +
    geom_point(size = 1) +
    geom_vline(
      aes(xintercept = 0),
      linetype = 2,
      colour = "black"
    ) +
    scale_colour_manual(values = gp_colours_full) +
    facet_wrap(
      ~ treat,
      nrow = 5,
      labeller = labeller(treat = label_intensity)
    ) +
    labs(x = "Time (week)", y = "Mean abundance (log)", colour = "Functional\ngroup") +
    estar:::theme_estar()
)

## ----echo=FALSE, message = FALSE, fig.height = 6, fig.width = 8, fig.cap="Organism count (mean +- sd in grey) of the five functional groups over the course of ecotoxicological experiment. The vertical line identifies the time when the insecticide was applied (at week 0), and facets represent the different concentrations (micro g/L)."----
(
  aquacomm_resps.rawplot <- aquacomm_fgps |>
    tidyr::pivot_longer(
      cols = c(herb, detr_herb, carn, omni, detr),
      names_to = "gp",
      values_to = "abund"
    ) |>
    ## summarize to plot
    dplyr::group_by(time, treat, gp) |>
    dplyr::summarize(abund.mean = mean(abund), abund.sd = sd(abund)) |>
    dplyr::ungroup() |>
    dplyr::mutate(
      gp = forcats::fct_recode(
        factor(gp),
        Herbivore = "herb",
        Detr_Herbivore = "detr_herb",
        Carnivore = "carn",
        Omnivore = "omni",
        Detritivore = "detr"
      )
    ) |>
    ggplot(aes(
      x = time,
      y = abund.mean,
      group = gp,
      colour = gp
    )) +
    geom_ribbon(
      aes(ymin = abund.mean - abund.sd, ymax = abund.mean + abund.sd),
      col = "gainsboro",
      fill = "gainsboro"
    ) +
    geom_line() +
    geom_point(size = 1) +
    geom_vline(
      aes(xintercept = 0),
      linetype = 2,
      colour = "black"
    ) +
    scale_colour_manual(values = gp_colours_full) +
    facet_wrap(
      ~ treat,
      nrow = 5,
      labeller = labeller(treat = label_intensity)
    ) +
    labs(x = "Time (week)", y = "Mean abundance", colour = "Functional\ngroup") +
    estar:::theme_estar()
)

## -----------------------------------------------------------------------------
invariability(
  type = "functional",
  mode = "lm_res",
  response = "lrr",
  metric_tf = c(1, max(aquacomm_resps$time)),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)

## -----------------------------------------------------------------------------
invariability(
  type = "functional",
  mode = "lm_res",
  metric_tf = c(1, max(aquacomm_resps$time)),
  response = "v",
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)

## -----------------------------------------------------------------------------
invariability(
  type = "functional",
  response = "lrr",
  mode = "cv",
  metric_tf = c(1, max(aquacomm_resps$time)),
  vd_i = "statvar_db",
  td_i = "time",
  b_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  d_data = aquacomm_resps
)

## -----------------------------------------------------------------------------
invariability(
  type = "functional",
  response = "v",
  mode = "cv",
  metric_tf = c(1, max(aquacomm_resps$time)),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)

## ----invar_schem, echo=FALSE, fig.height = 4, fig.cap="Figure 2: Schematic representation of how functional invariability is calculated by `estar`: $lrr$ and $diff$ refer to the log response ratio and the difference between the state variable in the disturbed system compared to the baseline time-series, $t_d$ is the time of disturbance. a) and b) illustrate the calculation of invariability as the coefficient of variation of the log response ratio and the difference, respectively; c and d) illustrate invariability calculated as the standard deviation of the residuals of the fitted linear model (blue solid line) of the log response ration and difference, respectively, as functions of time.", out.width = '500px'----
knitr::include_graphics("figures/invariability.png")

## -----------------------------------------------------------------------------
resistance(
  type = "functional",
  b = "input",
  res_mode = "lrr",
  res_time = "defined",
  res_t = 1,
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)

## -----------------------------------------------------------------------------
resistance(
  type = "functional",
  b = "input",
  res_mode = "lrr",
  res_time = "max",
  res_tf = c(0.14, 28),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)

## -----------------------------------------------------------------------------
resistance(
  type = "functional",
  b = "input",
  res_mode = "diff",
  res_time = "defined",
  res_t = 1,
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)

## -----------------------------------------------------------------------------
resistance(
  type = "functional",
  b = "input",
  res_mode = "diff",
  res_time = "max",
  res_tf = c(0.14, 28),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)

## -----------------------------------------------------------------------------
resistance(
  type = "functional",
  b = "d",
  b_tf = c(-4, 0.14),
  res_time = "max",
  res_mode = "diff",
  res_tf =c(0.14, 28),
  vd_i = "statvar_bl",
  td_i = "time",
  d_data = aquacomm_resps
)

## -----------------------------------------------------------------------------
resistance(
  type = "functional",
  b = "d",
  b_tf = c(-4, 0.14),
  res_mode = "lrr",
  res_time = "defined",
  res_t = 1,
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps
)

## ----resistance_schem, fig.height = 4, echo=FALSE, fig.cap="Figure 3: Schematic representation of how functional rate of recovery is calculated by `estar`: $lrr$ and $diff$ refer, respectively, to the log response ratio and to the difference between the state variable in the disturbed system in relation to the baseline value(s), $v$ to 'state variable', $t_{res}$, is the user-defined time step when resistance is calculated, and $t_{high}$, to the time step of the highest response in relation to the baseline. The four possibilities of calculating resistance (a-d) resistance result from calculating resistance from $lrr$ and $diff$ at time steps $t_{res}$ and $t_{high}$.", out.width = '500px'----
knitr::include_graphics("figures/resistance.png")

## -----------------------------------------------------------------------------
recovery_extent(
  type = "functional",
  response = "lrr",
  b = "input",
  t_rec = 28,
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
) 

## -----------------------------------------------------------------------------
recovery_extent(
  type = "functional",
  response = "lrr",
  b = "d",
  summ_mode = "mean",
  b_tf = c(-4, 0),
  t_rec = 28,
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps
)

## -----------------------------------------------------------------------------
recovery_extent(
  type = "functional",
  response = "diff",
  b = "input",
  t_rec = 28,
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)

## ----echo=FALSE, fig.height = 4, fig.cap="Figure 4: Schematic representation of how functional extent of recovery is calculated by `estar`. $t_{post}$ is the user-defined timestep when recovery is supposed to have happened. a) $lrr$ is the log response ratio of the state variable in the disturbed system compared to the baseline time-series and b) $diff$ to the difference between them.", out.width = '500px'----
knitr::include_graphics("figures/extent_recovery.png")

## -----------------------------------------------------------------------------
recovery_rate(
  type = "functional",
  response = "lrr",
  b = "input",
  metric_tf = c(0.14, 28),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)

## -----------------------------------------------------------------------------
recovery_rate(
  type = "functional",
  response = "v",
  b = "input",
  metric_tf = c(0.14, 28),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)

## ----echo=FALSE, fig.height = 4, fig.cap="Figure 5: Schematic representation of how functional rate of recovery is calculated by `estar`: a) $lrr$ is the log response ratio of the state variable in the disturbed system compared to the baseline time-series, b) $v$ is the state variable. The blue solid line identifies the linear model fitted to predict the response from time (along with its equation) and the double-headed arrow identifies the slope, which is the measure of rate of recovery", out.width = '500px'----
knitr::include_graphics("figures/rate_recovery.png")

## -----------------------------------------------------------------------------
persistence(
  type = "functional",
  b = "input",
  metric_tf = c(28, 56),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)

## -----------------------------------------------------------------------------
persistence(
  type = "functional",
  b = "d",
  b_tf = c(-4, 0.14),
  metric_tf = c(28, 56),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps
)

## ----echo=FALSE, out.width = '300px', fig.cap="Figure 6: Schematic representation of how functional persistence is calculated by `estar`: $v_b \\pm sd_b$ is the  standard deviation limit around the mean ($v_b$) of the baseline state variable. $t_a$ is the period for which persistence is calculated for, and $t_P$, the period during which the state variable remained inside the limits."----

knitr::include_graphics("figures/persistence.png")

## -----------------------------------------------------------------------------
oev(
  type = "functional",
  response = "lrr",
  metric_tf = c(0.14, 56),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)

## ----echo=FALSE, out.width = '300px', fig.cap="Figure 7: Schematic representation of the calculation of overall ecological etability in `estar`: $lrr$ is the log response ratio of the state variable in the disturbed system compared to the baseline time-series, $t_d$ is the time step when disturbance happened. The purple shaded area is the area under the curve calculated by the metric."----

knitr::include_graphics("figures/oev.png")

## -----------------------------------------------------------------------------
# Reformatting the macroinvertebrate community long-format data into 
# community composition data
comm_data <- aquacomm_fgps |>
    dplyr::group_by(time, treat) |>
    dplyr::filter(treat %in% c(0, 44)) |>
    dplyr::summarize_at(vars(herb, detr_herb, carn, omni, detr),
                        mean) |>
    dplyr::ungroup()

control_comm <- comm_data |>
        dplyr::filter(treat == 0) |>
        dplyr::select(-treat)
dist_comm <- comm_data |>
        dplyr::filter(treat == 44) |>
        dplyr::select(-treat)

## ----eval=FALSE---------------------------------------------------------------
# invariability(
#   type = "compositional",
#   metric_tf = c(0.14, 56),
#   comm_d = dist_comm,
#   comm_b = control_comm,
#   comm_t = "time")

## -----------------------------------------------------------------------------
resistance(type = "compositional",
          res_tf = c(0.14, 28),
          res_time = "max",
          comm_d = dist_comm,
          comm_b = control_comm,
          comm_t = "time")

## ----eval=FALSE---------------------------------------------------------------
# resistance(type = "compositional",
#           res_time = "defined",
#           res_t = 28,
#           comm_d = dist_comm,
#           comm_b = control_comm,
#           comm_t = "time")

## -----------------------------------------------------------------------------
recovery_rate(type = "compositional",
              metric_tf = c(0.14, 28),
              comm_d = dist_comm,
              comm_b = control_comm,
              comm_t = "time")

## -----------------------------------------------------------------------------
recovery_extent(type = "compositional",
          t_rec = 28,
          comm_d = dist_comm,
          comm_b = control_comm,
          comm_t = "time")

## -----------------------------------------------------------------------------
persistence(type = "compositional",
            b = "input",
            metric_tf = c(28, 56),
            comm_d = dist_comm,
            comm_b = control_comm,
            comm_t = "time",
            low_lim = 0.5,
            high_lim = 0.9)

## -----------------------------------------------------------------------------
oev(type = "compositional",
    metric_tf = c(0.14, 56),
    comm_d = dist_comm,
    comm_b = control_comm,
    comm_t = "time")

## ----echo=FALSE---------------------------------------------------------------
load("functional_performance.rda")

## -----------------------------------------------------------------------------
df_list$inv_benchmark %>%
  dplyr::select(-neval) %>% 
  kableExtra::kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

## -----------------------------------------------------------------------------
df_list$resis_benchmark %>%
  dplyr::select(-neval) %>% 
  kableExtra::kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

## -----------------------------------------------------------------------------
df_list$rate_benchmark %>%
  dplyr::select(-neval) %>% 
  kableExtra::kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

## -----------------------------------------------------------------------------
df_list$extent %>%
  dplyr::select(-neval) %>% 
  kableExtra::kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

## -----------------------------------------------------------------------------
df_list$persist_benchmark %>%
  dplyr::select(-neval) %>% 
  kableExtra::kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

## -----------------------------------------------------------------------------
df_list$oev_benchmark %>%
  dplyr::select(-neval) %>% 
  kableExtra::kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))