## ----knit-setup, include=FALSE, message=FALSE---------------------------------
knitr::opts_chunk$set(
  comment = "#>",
  warning=FALSE
  )
options(rmarkdown.html_vignette.check_title = FALSE) ## title and VignetteIndexEntry are slightly different and this avoids a warning message every time it's rendered

## ----message=FALSE------------------------------------------------------------
library(estar)
library(tidyr)
library(dplyr)
library(tibble)
library(ggplot2)
library(purrr)
library(viridis)
library(MARSS)
library(hesim)
library(kableExtra)
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 = " Figure 3: Organism mean count (log axis) of the five functional groups over the course of ecotoxicological experiment. The vertical dashed line identifies the time when chloryfiros insecticide was applied (at week 0), and facets represent the different concentrations (micro g/L)."----
(data_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())

## -----------------------------------------------------------------------------
Z_I5 <- matrix(list(0), 5, 5)
diag(Z_I5) <- 1

## ----collapse=TRUE------------------------------------------------------------
## calculate z-scores of abundances
data_zscores <- aquacomm_fgps %>%
  dplyr::filter(time == 0.14 | (time >= 4 && time <= 28)) %>%
  dplyr::mutate_at(vars(herb, detr_herb, carn, omni, detr),
                   ~if_else(. == 0, log(. + 0.001), log(.))) %>%
  ungroup() %>%
  group_by(., repli) %>%
  dplyr::mutate_at(vars(herb, detr_herb, carn, omni, detr),
                   ~MARSS::zscore(.)) %>%
  tidyr::pivot_longer(cols = c(herb, detr_herb, carn, omni, detr),
                      names_to = "fgp",
                      values_to = "abund_z")

## summarize abunds over replicates
data_zsumm <- data_zscores %>%
    dplyr::group_by(time, treat, fgp) %>%
    dplyr::summarize(abundz_mu = mean (abund_z),
                     abundz_sd = sd(abund_z)) %>%
    dplyr::ungroup()

## convert into time-series matrix
data_matrix <- data_zsumm %>%
    dplyr::select(time, treat, fgp, abundz_mu) %>%
    split(.$treat) %>%
    purrr::map(~ dplyr::select(., time, fgp, abundz_mu) %>%
        unique() %>%
        tidyr::pivot_wider(id_cols = fgp,
                           names_from = time, values_from = "abundz_mu", 
                           names_prefix = "time ") %>%
        tibble::column_to_rownames(var = "fgp") %>%
        as.matrix())

## ----collapse=TRUE------------------------------------------------------------
## 5 groups
### no observation error variance: 0 reps (summarized data), 5 gps
R_05 <- matrix(list(0), 5, 5)

data_marss <- data_matrix %>%
    purrr::map(~ MARSS(.,
                list(B = "unconstrained",
                     U = "zero", A = "zero",
                     Z = "identity", ## Z_I5 could also have been provided here 
                     Q = "diagonal and equal", 
                     R = R_05,
                     tinitx = 1),
                method = "BFGS"))
names(data_marss) <- paste0("Conc. = ", c("0", "0.1", "0.9", "6", "44" ), " micro g/L")

## -----------------------------------------------------------------------------
(data_Bs <- data_marss %>%
  purrr::map(~estar::extractB(., 
                states_names = c("Herbivores", "DetHerbivores", "Carnivores", "Omnivores", "Detrivores"))))

## -----------------------------------------------------------------------------
data_Bs %>%
  purrr::map(estar::reactivity)

## -----------------------------------------------------------------------------
data_Bs %>%
  purrr::map(estar::max_amp)

## -----------------------------------------------------------------------------
data_Bs %>%
  purrr::map(estar::init_resil)

## -----------------------------------------------------------------------------
data_Bs %>%
  purrr::map(estar::asympt_resil)

## -----------------------------------------------------------------------------
data_Bs %>%
  map(estar::stoch_var)

## -----------------------------------------------------------------------------
load("jacobian_performance.rda")

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

## ----sessionInfo--------------------------------------------------------------
Sys.time()
sessionInfo()