## ----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 the possible metrics of invariability available in the estar package. 'lrr' refers to the log response ratio of the state variable in the disturbed system compared to the baseline time-series, and 'v' to 'state variable' in the disturbed system. a) and b) demonstrate calculation of invariability based on the standard deviation of the residuals of the fitted linear model, whereby the linear model fit is shown by a blue solid line. c) and d) demonstrate calculation of invariability based on the coefficient of variation", out.width = '100%'----
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(1, 20),
  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(1, 20),
  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(1, 20),
  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 the possible metrics of resistance available in the estar package. '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', 'td+1', to the first time step after disturbance, and 'tlow', to the time step of the highest response in relation to the baseline.", out.width = '100%'----
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 = "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
)

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

## ----echo=FALSE, fig.height = 4, fig.cap="Figure 4: Schematic representation of the possible metrics of extent of recovery available in the estar package. 'lrr' refers to the log response ratio of the state variable in the disturbed system compared to the baseline time-series, 'diff' to the difference between them, and 'tpost', to a post-disturbance time step, specified by the user.", out.width = '100%'----
knitr::include_graphics("figures/extent_recovery.png")

## -----------------------------------------------------------------------------
recovery_rate(
  type = "functional",
  response = "v",
  b = "input",
  metric_tf = c(1, 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(1, 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 the possible metrics of rate of recovery available in the estar package. 'lrr' refers to the log response ratio of the state variable in the disturbed system compared to the baseline time-series, 'v' to '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 = '100%'----
knitr::include_graphics("figures/rate_recovery.png")

## -----------------------------------------------------------------------------
persistence(
  type = "functional",
  b = "input",
  metric_tf = c(28, 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
)

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

## ----echo=FALSE, out.width = '50%',vfig.cap="Figure 6: Schematic representation of the possible metrics of persistence available in the estar package. 'v' refers to 'state variable', 'lim_u' and 'lim_l' are the upper and lower limits defined as +1 and -1 standard deviation of the baseline state variable around the mean of the baseline state variable. 'ta' is the period for which persistence is calculated for, and 'tP', 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 = '50%',vfig.cap="Figure 7: Schematic representation of the calculation of Overall Ecological Stability in estar. 'lrr' refers to the log response ratio of the state variable in the disturbed system compared to the baseline time-series, 'td' is the time step where 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"))