## ----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 2: 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)."----
(aquacomm_fgps.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
aquacommz_allscen.ldf <- aquacomm_fgps %>%
    dplyr::filter(time >= 1 , time <= 28) %>%
    ungroup() %>%
    dplyr::mutate_at(vars(herb, detr_herb, carn, omni, detr), 
                     ~MARSS::zscore(.)) %>%
    dplyr::mutate(across(c(herb, detr_herb, carn, omni, detr),
                         ~dplyr::na_if(., 0))) %>%
    tidyr::pivot_longer(cols = c(herb, detr_herb, carn, omni, detr),
                        names_to = "fgp",
                        values_to = "abund_z")

## summarize abunds over replicates
aquacommz_allscen.summldf <- aquacommz_allscen.ldf %>%
    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
aquacommz_allscen.summmxls <- aquacommz_allscen.summldf %>%
    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)

aquacommz_allscen.marssls <- aquacommz_allscen.summmxls %>%
    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(aquacommz_allscen.marssls) <- paste0("Conc. = ", c("0", "0.1", "0.9", "6", "44" ), " micro g/L")

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

## -----------------------------------------------------------------------------
aquacomm.Bls %>%
  purrr::map(estar::reactivity)

## -----------------------------------------------------------------------------
aquacomm.Bls %>%
  purrr::map(estar::max_amp)

## -----------------------------------------------------------------------------
aquacomm.Bls %>%
  purrr::map(estar::init_resil)

## -----------------------------------------------------------------------------
aquacomm.Bls %>%
  purrr::map(estar::asympt_resil)

## -----------------------------------------------------------------------------
aquacomm.Bls %>%
  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()