## ----include = FALSE, echo = FALSE--------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup, echo = FALSE, include = FALSE-------------------------------------
library(ecorisk)


## ----indicator overview, echo = FALSE-----------------------------------------
ind_data <- data.frame(
  Trait = c("Cod", " Herring"), 
  General = c("Phytoplankton", "Seabirds"),
  Model = c("Cod spawning stock biomass", "Zooplankton mean size")
)

knitr::kable(ind_data, booktabs = TRUE) |>
  kableExtra::add_header_above(c("Expert scoring pathway" = 2, 
    "Modelling pathway" = 1)) 

## ----create tables, class.output="scroll-100"---------------------------------
exposure_scoring <- create_template_exposure(
  pressures = c("temperature", "salinity", "oxygen", "nutrient", "fishing"),
  n_components = 4, 
  mode_uncertainty = "component"
)
names(exposure_scoring)
# Rename exposure components
names(exposure_scoring)[2:9] <- c("magnitude", "frequency", "trend", "spatial",
                                  "uncertainty_magnitude", "uncertainty_frequency",
                                  "uncertainty_trend", "uncertainty_spatial")

sensitivity_scoring <- create_template_sensitivity(
  indicators = c("phytoplankton", "herring", "cod", "seabirds"),
  pressures = c("temperature", "salinity", "oxygen", "nutrient", "fishing"),
  type = c("direct", "direct_indirect"),
  n_sensitivity_traits = 5,
  mode_adaptive_capacity = "trait",
  mode_uncertainty = "trait"
)
names(sensitivity_scoring)

# Replaice the generic traits names ('...trait_1', '...trait_2') with
# the actual trait names
names(sensitivity_scoring) <- names(sensitivity_scoring) |> 
  stringr::str_replace("trait_1$", "feeding") |> 
  stringr::str_replace("trait_2$", "behaviour") |> 
  stringr::str_replace("trait_3$", "reproduction") |> 
  stringr::str_replace("trait_4$", "habitat") |> 
  stringr::str_replace("trait_5$", "general") 

## ----read expert scores-------------------------------------------------------
exposure_scoring <- ex_expert_exposure
head(exposure_scoring)
sensitivity_scoring <- ex_expert_sensitivity
head(sensitivity_scoring)

## ----read time series data----------------------------------------------------
ts_pressures <- pressure_ts_baltic
head(ts_pressures)
ts_indicators <- indicator_ts_baltic
head(ts_indicators)

## ----exposure and sensitivity expert pathway, echo = TRUE---------------------
exp_expert <- calc_exposure(
  pressures = exposure_scoring$pressure, 
  components = exposure_scoring[ ,2:5],
  uncertainty = exposure_scoring[ ,6:9],
  method = "mean"
)
head(exp_expert)

sens_ac_expert <- calc_sensitivity(
  indicators = sensitivity_scoring$indicator, 
  pressures = sensitivity_scoring$pressure,
  type = sensitivity_scoring$type,
  sensitivity_traits = sensitivity_scoring[ ,4:8],
  adaptive_capacities = sensitivity_scoring[ ,9:13],
  uncertainty_sens = sensitivity_scoring[ ,14:18],
  uncertainty_ac = sensitivity_scoring[ ,19:23], 
  method = "mean"
)
head(sens_ac_expert)

## ----model exposure-----------------------------------------------------------
exposure_model <- model_exposure(
  pressure_time_series = ts_pressures,
  base_years = c(start = 1984, end = 1993),
  current_years = c(start = 2007, end = 2016),
  trend = "return", 
  spatial = c(2, 2, 5, 5, 3, 3, 2, 2)
)
exposure_model

## ----model sensitivity--------------------------------------------------------
sens_ac_model <- model_sensitivity(
  indicator_time_series = ts_indicators, 
  pressure_time_series = ts_pressures,
  current_years = c(start = 2010, end = 2016)
)
sens_ac_model

## ----model diagnostics sensitivity--------------------------------------------
plot_diagnostic_sensitivity(
   indicator_time_series = ts_indicators[, c(1:2)],
   pressure_time_series = ts_pressures[, c(1:2)]
)

## ----vulnerability------------------------------------------------------------
vuln_experts <- vulnerability(
  exposure_results = exp_expert, 
  sensitivity_results = sens_ac_expert,
  method_vulnerability = "mean",
  method_uncertainty = "mean"
)

vuln_model <- vulnerability(
  exposure_results = exposure_model, 
  sensitivity_results = sens_ac_model,
  method_vulnerability = "mean",
  method_uncertainty = "mean"
)

head(vuln_experts)
head(vuln_model)

## ----status experts-----------------------------------------------------------
status_experts <- data.frame(
  indicator = c("phytoplankton", "herring", "cod", "seabirds"), 
  status = c("good", "undesired", "undesired", "good"),
  score = c(1, -1, -1, 1)
)
status_experts

## ----status model-------------------------------------------------------------
status_model <- status(
  indicator_time_series = ts_indicators,
  base_years = c(start = 1984, end = 1993),
  current_years = c(start = 2012, end = 2016),
  sign = "-", 
  condition = "<"
)
status_model

## ----risk---------------------------------------------------------------------
risk_expert <- risk(
  vulnerability = vuln_experts, 
  status = status_experts
)
risk_model <- risk(
  vulnerability = vuln_model, 
  status = status_model
)

head(risk_expert)
head(risk_model)

## ----rename and select pressure variables model pathway-----------------------
risk_model <- risk_model[c(1, 3, 5, 7, 8, 9, 12, 14:16), ]
risk_model$pressure <- c(
  "nutrient", "temperature", "salinity", "oxygen", "fishing",   # for zooplankton
  "nutrient", "temperature", "salinity", "oxygen", "fishing")   # for cod

## ----aggregate pathways and risks---------------------------------------------
risks <- rbind(risk_expert, risk_model)

aggregated_risk <- aggregate_risk(
  risk_results = risks, 
  method = "mean"
)

aggregated_risk$multi_indicator_risk
aggregated_risk$multi_pressure_risk

## ----plot radar, fig.width=6, fig.height=6------------------------------------
p_radar <- plot_radar(
  risk_scores = risks,
  aggregated_scores = aggregated_risk,
  type = "direct_indirect", 
  pathway = "combined"
)

p_radar[[1]]
p_radar[[2]]
p_radar[[3]]
p_radar[[4]]
p_radar[[5]]
p_radar[[6]]

## ----correlation plot, fig.width=6, fig.height=6------------------------------
temp <- risks[c(26:30, 45:49), c(1, 2, 7)]
temp <- temp |>
  tidyr::pivot_wider(names_from = indicator, values_from = risk)

ggplot2::ggplot(dat = temp, 
    ggplot2::aes(x = cod, y = eastern_baltic_cod, colour = pressure)) +
  ggplot2::geom_point() +
  ggplot2::geom_abline(slope = 1, intercept = 0) +
  ggplot2::xlim(-10, 10) + 
  ggplot2::ylim(-10,10) +
  ggplot2::labs(x = "Expert-based pathway", y = "Modelling-based pathway") +
  ggplot2::theme_minimal()

## ----plot heatmap, fig.width = 8, fig.height=9, out.height="25%", out.width="%"----
p_heat <- plot_heatmap(
  risk_scores = risks,
  aggregated_scores = aggregated_risk,
  uncertainty = TRUE
)

p_heat[[1]]
p_heat[[2]]