oldHTM.Rmd

---
title: "Example #1: Calibration Protocol - Time-Averaged"
author: Julia L. Blanchard
date: July 2020
place: Hobart, Australia
output:
  html_document: default
  # pdf_document: default
always_allow_html: true
#runtime: shiny
---

# Calibrating a multi-species model to time-averaged species' catches

In this example we will explore how to calibrate a size spectrum model to data using the "mizer" R package. 

Recall, there are three different kinds of size spectrum models in mizer, of increasing complexity: 

1) community model: purely size-based and representative of a single but  "average" species across the whole community.

2) trait-based model, which disaggregates the size spectrum into different groups with different life-histories, through differences in each "species" asymptotic which determines other life-history parameters such as the size at maturity (Hartvig et al. 2011, Andersen & Pedersen, 2010).

3) multispecies model - which has the same equations and parameters as the trait-based model but is parameterised to represent multiple species in a real system, where each species can have many differing species-specific traits (Blanchard et al. 2014). 

Here we focus on multispecies size spectrum models. In practice, these models have been parametrised in a few different ways depending on data availability for a system or research questions.

Some studies have focused on  many species-specific values, for example where each species have different values of life-history, size-selective feeding trait parameters (e.g. $\beta$ and $\sigma$), and details of species interactions (Blanchard et al. 2014, Reum et al. 208) to better capture the dynamics of marine food webs. 

Others, such as Jacobsen et al. (2014, 2016), have represented variation in only a couple of the most important life history parameters for each species - asymptotic size (which links to other parameters such as maturation size), growth, vulnerability and recruitment parameters ($R_{max}$, $eRepro$) to broadly capture fished communities or carry out across ecosystem comparisons. 

Once you have parametrised the multispecies model for your system (section 1), you may find that species do not coexist or the biomass or catches are very different from your observations. After the model is parameterised and assessed for basic principles and coexistence (section 3), further calibration to observational data is used to ensure the relative abundance of each species reflects the system (section 4), at least during a stable period, which is time-averaged. 

The background resource parameters and the recruitment parameters, $R{max}$ (maximum recruitment) and $erepro$ (reproductive efficiency) greatly affect the biomasses of species in your system. 

The recruitment parameters capture those density dependent effects that are not explicitly modelled. Here, we use a Beverton-Holt type stock recruitment relationship to represent these effects.

As a starting point, we will estimate these parameters as a means of fitting the modelled species catches to the observed catches. This could similarly be carried out with biomasses. Other model detailed approaches also exist, see the main paper, but this approach has been used to get models in the right "ball-park", which can them be further evaluated using diagnostics (example X) and fitted to time series data (example XX).




### A Simple Protocol for Multispecies Model Calibration

We will adapt the "recipe" for calibration in Jacobsen et al 2014 (see supp. mat.) and Blanchard et al (2014), into the following steps:

0. Run the model with the chosen species-specific parameters. This will relate some of the missing parameters to $w_{inf}$ ($h$ and $\gamma$ - explain through simple example of how the model works?). $R_{max}$ (see example that explains $R_{max}$?) could also be automatically calculated based on equilbrium assumptions (Andersen et al. 2016) but by default it is "$Inf$", which means there is no density dependence associated with spawner-recruit dynamics (RF: default is the 2016 method at the moment but setting up at $Inf$ to start with something that doesn't coexist below).

1. Obtain the time-averaged data (e.g. catches or biomasses for each species) and the time-averaged fishing mortalty inputs (e.g. from stock assessments). Typically this should be over a stable part of the time series for your system.

2. Calibrate the maximum recruitment, $R_{max}$ by hand by penalising the reproductive capacity of the largest species.  Check that the physiological recruitment, $RDI$, is much higher than the realised recruitment, $RDD$.

3. Calibrate the maximum recruitment, $R_{max}$ with optimisation, by minimising the error between observed and estimated catches. 

4. Refine the calibration of the recruitment via $eRepro$ by hand and check that the species' $F_{msy}$ are reasonable.

5. Calibrate growth by changing the carrying capacity of the background resource spectrum ($\kappa$).

6. Calibrate the food competition


## Step 0. Run the model with the chosen species-specific parameters.

In this section you will:

- obtain or create a dataframe of species-specific parameters

- run the dataframe through the `mizer` package and examine the model output


A species-specific dataframe is already stored in `mizer`, which contains the North Sea Model Parameters (RF however is probably not updated so using the .csv in the repository)


```{r step 0 - initial parameters, message= F }
# if user has not installed the requird packages
# install.packages("tidyverse")
# install.packages("plotly")
# devtools::install_github("sizespectrum/mizer")
#devtools::install_github("sizespectrum/mizerExperimental")
library(mizerExperimental) # for projectToSteady()
library(mizer)
library(tidyverse)
library(plotly)
library(tictoc)
library(shiny)
library(shinyWidgets)
library(parallel)
library(optimParallel)
library(ggrepel)
# loading North Sea data
nsParams <- read.csv("data/nsparams.csv")
# This data frame already has Rmax values, let's remove them to calibrate them again later
nsParams[,"r_max"] <- Inf
# ordering by asymptotic size for color gradient
nsParams <- nsParams[order(nsParams$w_inf),] 
# need to order the interaction matrix as well | it's so ugly please help
inter <- select(as.data.frame(inter),nsParams$species)
inter <- as.matrix(select(as.data.frame(t(inter)),nsParams$species))
# If you want to make it less multi-species and more trait-based model
# nsParams[,"beta"] <-100
# nsParams[,"sigma"] <-1.5
```


```{r step 0 - loading plot function, include = F}
require(cowplot)
require(gridExtra)
g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}
# arrow to show asymptotic size
plotSummary <- function (x, y, power = 2, wlim = c(.1,NA), short = F, ...) 
{
  xlim = c(NA,10^log10(max(x@params@species_params$w_inf)))
  font_size = 7
  
  # need to display the legend at the bottom and only p1 has the background so using that one
  
  p1 <- plotSpectra(x, power = power, wlim = wlim, ...)
  p1 <- p1  + scale_x_continuous(limits = xlim, trans = "log10", name = "Individual size [g]") +
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          text = element_text(size=font_size),
          panel.background = element_blank(),
          panel.grid.minor = element_line(color = "gray"),
          legend.position = "right", legend.key = element_rect(fill = "white"))
  # guides(color = guide_legend(nrow=2))
  
  mylegend<-g_legend(p1) # save the legend
  p1 <- p1 + theme(legend.position = "none") # now remove it from the plot itself
  
  p7 <- plotBiomass(x)
  p7 <- p7 + theme(legend.position = "none",
                   text = element_text(size=font_size),
                   panel.background = element_blank(),
                   panel.grid.minor = element_line(color = "gray"),)
  # theme_bw()
  
    if(short)
{
      p1 <- p1 +theme(axis.title.x=element_text(),
                    axis.text.x=element_text(),
                    axis.ticks.x=element_line())
      leftCol <- plot_grid(p1,p7,
                           ncol = 1, align = "v", axis = "l")
      p10 <- plot_grid(leftCol, mylegend,
                       rel_widths = c(6,1),
                       ncol = 2)
  
} else
{
  
  
  p2 <- plotFeedingLevel(x, include_critical = F, ...)
  p2 <- p2 + scale_x_continuous(limits = xlim, trans = "log10") +
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          text = element_text(size=font_size),
          legend.position = "none")
  p3 <- plotPredMort(x, ...)
  p3 <- p3 + scale_x_continuous(limits = xlim, trans = "log10") +
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          text = element_text(size=font_size),
          legend.position = "none")
  p4 <- plotFMort(x, ...)
  p4 <- p4 + scale_x_continuous(limits = xlim, trans = "log10", name = "Individual size [g]") +
    theme(legend.position = "none",
          text = element_text(size=font_size),)
  
  
  
  
# yeild and ssb |
  
  bm <- getBiomassFrame(x)
  plot_dat <- filter(bm, Year == max(unique(bm$Year)))
  # plot_dat$w_inf <- x@params@species_params$w_inf
  yieldDat <- getYield(x)
  plot_dat$yield <- yieldDat[dim(yieldDat)[1],]
  plot_dat$Year <- NULL
  plot_dat <- melt(plot_dat,"Species")
  p5 <- ggplot(plot_dat) +
    geom_bar(aes(x = Species,y = value, fill = Species, alpha = variable), stat = "identity", position = position_dodge()) +
    coord_cartesian(ylim = c(0.5*min(plot_dat$value),NA)) +
    # geom_text(aes(x = Species, y = value, label = Species), check_overlap = T)+
    
    # geom_point(aes(x = w_inf, y = Biomass, color = species)) +
    # geom_point(aes(x = w_inf, y = yield, color = species), shape = "+", size = 5) +
    # scale_x_continuous(name = "SSB and Yield") +
    scale_y_continuous(trans = "log10", name = "SSB and Yield") + #, limits = c(0.5*min(plot_dat$value),NA)) +
    scale_fill_manual(name = "Species", values = x@params@linecolour) +
    scale_alpha_manual(name = "Stat", values = c(1, 0.5), labels = c("SSB","Yield")) +
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          text = element_text(size=font_size),
          legend.position = "bottom", legend.key = element_rect(fill = "white"))
  
  mylegend<-g_legend(p5) # save the legend
  p5 <- p5 + theme(legend.position = "none")
  
  # try yield divided by total biomass, we want to see the difference
  
  # r0
  
  plot_dat <- as.data.frame(getRDI(x@params)/getRDD(x@params))
  plot_dat$species <- factor(rownames(plot_dat),x@params@species_params$species)
  colnames(plot_dat)[1] <- "ratio"
  plot_dat$w_inf <- as.numeric(x@params@species_params$w_inf)
  
  # trying to have bars at their w_inf but on a continuous scale
  plot_dat$label <- plot_dat$species
  plot_dat2 <- plot_dat
  plot_dat2$ratio <- 0
  plot_dat2$label <- NA
  plot_dat <- rbind(plot_dat,plot_dat2)
  
  
  p6 <- ggplot(plot_dat) +
    geom_line(aes(x = w_inf, y = ratio, color = species), size = 15, alpha = .8) +
    geom_text(aes(x = w_inf, y = ratio, label = label),position = position_stack(vjust = 0.5), angle = 30, size = 3)+
    scale_color_manual(name = "Species", values = x@params@linecolour) +
    scale_y_continuous(name = "RDI/RDD") +
    scale_x_continuous(name = "Asymptotic size (g)", trans = "log10") +
    theme(legend.position = "none",
          text = element_text(size=font_size),) 
  
  p10 <- plot_grid(p1,p2,p3,p4, p5, p6, p7, mylegend, byrow = F,  
                   # rel_heights = c(1,1,1,2),
                   rel_widths = c(2,1),
            nrow = 4, align = "v", axis = "l")
}
 # p <- grid.arrange(p10,mylegend, nrow=2,heights=c(9.5,0.5))
return(p10)
}
# hopefully all of this will go on sizespectrum/mizer, in the meantime
plotPredObsYield <-function(sim, dat, returnData = FALSE){
## check obs vs. predicted yield
plot_dat <-melt(getYield(sim)[100,]/1e6)
plot_dat$obs <- log10(dat)
plot_dat$value <- log10(plot_dat$value)
plot_dat$Species <-row.names(plot_dat)
w_inf <- log10(sim@params@species_params$w_inf)
names(w_inf) <- sim@params@species_params$species
# window size
winLim <- c(min(plot_dat$obs,plot_dat$value), max(plot_dat$obs,plot_dat$value))
p <- ggplot(plot_dat) + # plot predicted and observed yields
  geom_point(aes(x = value, y = obs, color = Species, size = Species)) +
  geom_text_repel(aes(x = value, y = obs, label = Species), hjust = 0, nudge_x = 0.05)+
  scale_size_manual(values = w_inf) +
  scale_color_manual(values = sim@params@linecolour) +
  geom_abline(color = "black", slope = 1, intercept = 0, linetype = "dashed", alpha = .5) + 
  scale_x_continuous(name = "log10 Predicted Yield", limits = winLim) + 
  scale_y_continuous(name = "log10 Observed Yield", limits = winLim) +
  theme(legend.position = "none", legend.key = element_rect(fill = "white"),
        panel.background = element_blank(), panel.grid.minor = element_line(color = "gray"))
if(returnData) return(plot_dat) else return(p)
} 
plotDiet2 <- function (sim, species = NULL, xlim = c(1,NA), returnData = F) 
{
  params <- sim@params
    # if (is.integer(species)) {
    #     species <- params@species_params$species[species]
    # }
  
    # diet <- getDiet(params)[params@species_params$species == 
    #     species, , ]
    # prey <- dimnames(diet)$prey
    # prey <- factor(prey, levels = rev(prey))
    # plot_dat <- data.frame(Proportion = c(diet), w = params@w, 
    #     Prey = rep(prey, each = length(params@w)))
    # plot_dat <- plot_dat[plot_dat$Proportion > 0, ]
    # 
    # ggplot(plot_dat) + geom_area(aes(x = w, y = Proportion, fill = Prey)) + 
    #     scale_x_log10(limits = xlim) + labs(x = "Size [g]") + 
    #   scale_fill_manual(values = sim@params@linecolour) +
    #   ggtitle(species)
    
    
    diet <- getDiet(params)
    plot_dat <- melt(diet)
    plot_dat <- plot_dat[plot_dat$value > 0, ]
    colnames(plot_dat) <- c("Predator", "size", "Prey", "Proportion")
    
    if(is.null(species)) p <- ggplot(plot_dat) + facet_wrap(.~Predator, scales = "free") else p <- ggplot(filter(plot_dat, Predator == species))
    
  p <- p +
      geom_area(aes(x = size, y = Proportion, fill = Prey))+
       scale_x_continuous(limits = c(1,NA), name = "Size [g]", trans = "log10") + 
      scale_fill_manual(values = sim@params@linecolour)+
      theme(legend.position = "right", legend.key = element_rect(fill = "white"),
        panel.background = element_blank(), panel.grid.minor = element_line(color = "gray"),
        strip.background = element_blank())
  
  if(returnData) return(plot_dat) else return(p)
    
}
# Try facets
plotFmsy <- function(params, effortRes = 20, returnData = F, speciesData = NULL)
{
# make one gear per species so we can bary the effort per species
gear <- gear_params(params)
gear$gear <- params@species_params$species
gear_params(params) <- gear
catchability <- params@species_params$catchability
xlim <- 1.5 # maximum effort* catchability / xaxis limit
# we want to vary effort value so we get a scale from 0 to 1 of effort * catchability per species
# the "species" arg allows to run the function for only one species, which should be faster but it means "species" must also contain the result of every other species (so it's a two object list)
if(!is.null(speciesData))
{
  speciesName <- speciesData[[1]] # which species are we changing?
  plot_dat <- speciesData[[2]] # plot_dat of all species
  plot_dat <- filter(plot_dat, species!= speciesName) # remove previous result of the concerned species
  iSpecies <- which(params@species_params$species == speciesName)
  counter = 0 # sim counter
  # determine effort range
  effortMax <- round(xlim/catchability[iSpecies],1)+.1
  SpDat <- NULL
  
  effortSeq <- exp(seq(0,log(effortMax+1), length.out =  effortRes)) -1 
effortSeq <- effortSeq[effortSeq<effortMax] # creating an exponentially increasing effort sequence
  for(iEffort in effortSeq)
  {
    effort_vec <- rep(1,dim(params@species_params)[1]) # all effort set to one
    effort_vec[iSpecies] <- iEffort # except that one which varies
    if(!counter )
    {
       tempSim <- project(params, effort = effort_vec, t_max = 20)
       counter <- 1
    } else  tempSim <- project(params, effort = effort_vec, t_max = 10, initial_n = tempSim@n[dim(tempSim@n)[1],,], 
                       initial_npp = tempSim@n_pp[dim(tempSim@n_pp)[1],])
   #catch
  yieldDat <- getYield(tempSim)
  SpDat <- rbind(SpDat,c(yieldDat[dim(yieldDat)[1],iSpecies],iEffort))
  }
  SpDat <- as.data.frame(SpDat)
  SpDat$species <- params@species_params$species[iSpecies]
  SpDat$V2 <- SpDat$V2*catchability[iSpecies] # so V2 is effort * catchability
  colnames(SpDat) <- c("yield","effort","species")
plot_dat <- rbind(plot_dat,SpDat)
  
  
  
} else {
plot_dat <- NULL
for(iSpecies in 1:dim(params@species_params)[1])
{
  counter = 0 # sim counter
  # determine effort range
  effortMax <- round(xlim/catchability[iSpecies],1)+.1
  SpDat <- NULL
  
  effortSeq <- exp(seq(0,log(effortMax+1), length.out =  effortRes)) -1 # every .1 takes 2 min to run, evry .2 takes 1 min but lesser resolution
effortSeq <- effortSeq[effortSeq<effortMax] # creating an exponentially increasing effort sequence
  for(iEffort in effortSeq)
  {
    effort_vec <- rep(1,dim(params@species_params)[1]) # all effort set to one
    effort_vec[iSpecies] <- iEffort # except that one which varies
    if(!counter )
    {
       tempSim <- project(params, effort = effort_vec, t_max = 20)
       counter <- 1
    } else  tempSim <- project(params, effort = effort_vec, t_max = 10, initial_n = tempSim@n[dim(tempSim@n)[1],,], 
                       initial_npp = tempSim@n_pp[dim(tempSim@n_pp)[1],])
   #catch
  yieldDat <- getYield(tempSim)
  SpDat <- rbind(SpDat,c(yieldDat[dim(yieldDat)[1],iSpecies],iEffort))
  }
  SpDat <- as.data.frame(SpDat)
  SpDat$species <- params@species_params$species[iSpecies]
  SpDat$V2 <- SpDat$V2*catchability[iSpecies] # so V2 is effort * catchability
  colnames(SpDat) <- c("yield","effort","species")
plot_dat <- rbind(plot_dat,SpDat)
    
}
}
plot_dat$species <- factor(plot_dat$species, levels = params@species_params$species)
# colnames(plot_dat) <- c("yield","effort","species")
if(!is.null(speciesData)) p <- ggplot(filter(plot_dat, species == speciesName)) else p <- ggplot(plot_dat)
p <- p + geom_line(aes(x = effort , y = yield, color = species))+
  facet_wrap(species~., scales = "free") +
  scale_x_continuous(limits= c(0,xlim),name = "fishing mortality rate")+#, limits = c(1e10,NA))+
  scale_y_continuous(trans = "log10") +
  scale_color_manual(name = "Species", values = params@linecolour) +
    theme(legend.position = "none", legend.key = element_rect(fill = "white"),
           panel.background = element_blank(), panel.grid.minor = element_line(color = "gray"),
          strip.background = element_blank())
if(returnData) return(plot_dat) else return(p)
}
plotGrowthCurves2 <- function (object, 
                               species = NULL, 
                               max_age = 20, 
                               percentage = FALSE, 
                               species_panel = FALSE, 
                               highlight = NULL,
                               returnData = F) 
{
  if (is(object, "MizerSim")) {
    params <- object@params
    t <- dim(object@n)[1]
    params@initial_n[] <- object@n[t, , ]
    params@initial_n_pp <- object@n_pp[t, ]
  }
  else if (is(object, "MizerParams")) {
    params <- validParams(object)
  }
  species <- valid_species_arg(params, species)
  ws <- getGrowthCurves(params, species, max_age, percentage)
  plot_dat <- reshape2::melt(ws)
  plot_dat$Species <- factor(plot_dat$Species, params@species_params$species)
  plot_dat$legend <- "model"
  if (all(c("a", "b", "k_vb") %in% names(params@species_params))) {
    if ("t0" %in% names(params@species_params)) {
      t0 <- params@species_params$t0
    }
    else {
      t0 <- 0
    }
    VBdf <- data.frame(species = params@species_params$species, 
                       w_inf = params@species_params$w_inf, a = params@species_params$a, 
                       b = params@species_params$b, k_vb = params@species_params$k_vb, 
                       t0 = t0)
    VBdf$L_inf <- (VBdf$w_inf/VBdf$a)^(1/VBdf$b)
    plot_dat2 <- plot_dat
    plot_dat2$value <- apply(plot_dat, 1, function(x) {
      sel <- VBdf$species == x[1]
      length <- VBdf$L_inf[sel] * (1 - exp(-VBdf$k_vb[sel] * 
                                             (as.numeric(x[2]) - VBdf$t0[sel])))
      VBdf$a[sel] * length^VBdf$b[sel]
    })
    plot_dat2$legend <- "von Bertalanffy"
    plot_dat <- rbind(plot_dat, plot_dat2)
  }
  p <- ggplot(filter(plot_dat, legend == "model")) + geom_line(aes(x = Age, 
                                                                   y = value, colour = Species, linetype = Species, size = Species))
  y_label <- if (percentage) "Percent of maximum size" else "Size [g]"
  linesize <- rep(0.8, length(params@linetype))
  names(linesize) <- names(params@linetype)
  linesize[highlight] <- 1.6
  p <- p + scale_x_continuous(name = "Age [Years]") + scale_y_continuous(name = y_label) + 
    scale_colour_manual(values = params@linecolour) + scale_linetype_manual(values = params@linetype) + 
    scale_size_manual(values = linesize)
  if (!percentage) {
    if (length(species) == 1) {
      idx <- which(params@species_params$species == species)
      w_inf <- params@species_params$w_inf[idx]
      p <- p + geom_hline(yintercept = w_inf, colour = "grey") + 
        annotate("text", 0, w_inf, vjust = -1, label = "Maximum")
      w_mat <- params@species_params$w_mat[idx]
      p <- p + geom_hline(yintercept = w_mat, linetype = "dashed", 
                          colour = "grey") + annotate("text", 0, w_mat, 
                                                      vjust = -1, label = "Maturity")
      if ("von Bertalanffy" %in% plot_dat$legend) 
        p <- p + geom_line(data = filter(plot_dat, legend == 
                                           "von Bertalanffy"), aes(x = Age, y = value))
    }
    else if (species_panel) {
      p <- ggplot(plot_dat) + 
        geom_line(aes(x = Age, y = value, colour = legend)) + 
        scale_x_continuous(name = "Age [years]") +
        scale_y_continuous(name = "Size [g]") +
        facet_wrap(.~Species, scales = "free") +
        geom_hline(aes(yintercept = w_mat),
                   data = tibble(Species = as.factor(object@params@species_params$species[]),
                                 w_mat = object@params@species_params$w_mat[]),
                   linetype = "dashed", colour = "grey") +
        geom_hline(aes(yintercept = w_inf),
                   data = tibble(Species = as.factor(object@params@species_params$species[]),
                                 w_inf = object@params@species_params$w_inf[]),
                   linetype = "solid", colour = "grey") +
              theme(panel.background = element_blank(), panel.grid.minor = element_line(color = "gray"),
              strip.background = element_blank(), legend.key = element_blank())+
        scale_color_discrete(name = "Growth", labels = c("Modelled","von Bertalanffy"))
    }
  }
  if(returnData) return(plot_dat) else return(p)
}
```

```{r step 0 - parameters}
params_uncalibrated <- newMultispeciesParams(nsParams, inter, kappa = 1e11, max_w=1e6) # inter comes with loading "mizer"
params_uncalibrated@species_params$erepro <- .01 # default is 1
# note the volume of this model is set to the reflect the entire volume of the North Sea - hence the very large kappa value. This is system specific and you may wnat to work with per m^3 as in the defaults.
#  Add other params for info
#  param$Volumecubicmetres=5.5e13    #unit of volume. Here total volume of North sea is used (Andersen & Ursin 1977)
# have a look at species parameters that have been calculated
# params_uncalibrated@species_params
# alternative params without redundant parameters to reduce the size of the dataframe on the screeen
params_uncalibrated@species_params[,-which(colnames(params_uncalibrated@species_params) %in% 
                                             c("sel_func","gear","interaction_resource","pred_kernel_type","m","alpha","n","p","q","w_min"))]
```

w_inf: asymptotic size

w_mat: maturation size (determines when 50% of the population has matured / not sure!)

beta: preferred predator/prey mass ratio

sigma: width of the feeding kernel

R_max: Beverton-Holt density dependence parameter

k_vb: von Bertalanffy growth parameter

l25: length at ...

l50: length at ...

a: coefficient for age to size conversion

b: constant for age to size conversion

catchability: fisheries efficiency

h: maximum intake rate

k: metabolism constant

ks: metabolism coefficient

z0: background mortality coefficient

gamma: search volume (obtained from beta and sigma)

w_mat25: weight at which 25% of individuals are mature

erepro: coefficent that weights reproductive output


```{r step 0 - first run,message=F, warning = F}
# params_uncalibrated@species_params <- params_uncalibrated@species_params[order(params_uncalibrated@species_params$w_inf),] # ordering by asymptotic size for color gradient
# lets' change the plotting colours, by far the hardest and trickiest part
# looks good but hard to distinguish some species
library(viridis)
params_uncalibrated@linecolour[1:12] <-viridis(dim(params_uncalibrated@species_params)[1])
# easier to read plots but color meaningless
# library(pals)
# params_uncalibrated@linecolour[1:12] <-glasbey(dim(params_uncalibrated@species_params)[1])
# Gradient over asymptotic size but might get difficult to distinguish species past 9
# colfunc <- colorRampPalette(c("firebrick3","darkturquoise", "orange"))
# colGrad <- colfunc(dim(params_uncalibrated@species_params)[1])
#  #names(colGrad) <- params_uncalibrated@species_params$species[order(params_uncalibrated@species_params$w_inf)] # gradient is ordered per winf
# 
# params_uncalibrated@linecolour[1:12] <- colGrad
# 
params_uncalibrated@linecolour["Resource"] <-"seagreen3"
# jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000","maroon4"))
# params_uncalibrated@linecolour[1:12] <-jet.colors(dim(params_uncalibrated@species_params)[1])
# write color routine for color gradient readable but related to winf (can use my usual one but limited to ~9 species for clarity
# maybe have different color depening on species number
# line width related to asymptotic size ?
# run with fishing
sim_uncalibrated <- project(params_uncalibrated, t_max = 100, effort = 1)
plotSummary(sim_uncalibrated, short = T)
```

The top panel shows the different species size spectrum at the last time step of the simulation while the bottom panel shows the abundance per species through time.
These plots show that species do not coexist and several go extinct. This is because there was no external density dependence ($R_{max}$ is set at $Inf$) and the largest species (Cod and Saithe) are out-competing the rest.

## Step 1. Obtain the time-averaged data (e.g. catches or biomasses for each species) and the time-averaged fishing mortalty inputs (e.g. from stock assessments).

In this section you will:

- Download fisheries data and process them in a format comparable to the model output

- Visualise the data


The following .csv are extracted from the ICES database using "data/getICESFishdata_param.R". Fishing data is averaged over 2014-2019 as it is a relatively stable period in catches.

```{r step 1 - loading data}
# fisheries mortality F
fMat <- read.csv("data/fmat.csv")
# fMatWeighted <- read.csv("data/fmatWeighted.csv") # Sandeel and Cod have multiple data base so average their F and weighting by SSB
fMatWeighted <- readRDS("data/FmatWeightedInterpolated.rds") # to get Gurnard data
# read in time-averaged  catches  
catchAvg <-read.csv("data/time-averaged-catches.csv") # only that one is used at the moment | catches are estimated from fMatW
# ssb
ssbAvg <- read.csv("data/time-averaged-SSB.csv")
```



```{r step 1 - visual representation of the data, warning=F}
# plot_dat <- reshape2::melt(fMatWeighted,"X")
# colnames(plot_dat) <- c("Time", "Species", "F")
# 
# ggplot(plot_dat)+
#   geom_line(aes(x = Time, y = F, color = Species))+
#   scale_y_continuous(name = "Fisheries mortality rate") +
#   scale_color_manual(values = params_uncalibrated@linecolour) +
#   theme(                   panel.background = element_blank(),
#                    panel.grid.minor = element_line(color = "gray"),
#                    legend.key = element_rect(fill = "white"))
# overlapping fisheries mortality with next plot
plot_dat <- apply(as.matrix(fMatWeighted[67:73,]),2,mean, na.rm = T)
plot_dat <- plot_dat[2:length(plot_dat)]
plot_dat <- melt(plot_dat)
plot_dat$Species <- factor(as.character(rownames(plot_dat)),levels = c(as.character(params_uncalibrated@species_params$species)))
plot_dat <- plot_dat[order(plot_dat$Species),]
plot_dat$w_inf <- params_uncalibrated@species_params$w_inf
plot_dat2 <- plot_dat
  
plot_dat <- data.frame(catchAvg,ssbAvg)
plot_dat$species.1 <- NULL
colnames(plot_dat) <- c("Species", "average catch", "average SSB")
plot_dat$Species <- factor(as.character(plot_dat$Species),levels = c(as.character(params_uncalibrated@species_params$species)))
plot_dat <- reshape2::melt(plot_dat,"Species")
plot_dat$w_inf <- rep(params_uncalibrated@species_params$w_inf,2)
  
  # ggplot(plot_dat) +
  #   geom_bar(aes(x = Species,y = value, fill = Species, alpha = variable), stat = "identity", position = position_dodge()) +
  #   geom_point(aes(x = Species))
  #   scale_y_continuous(trans = "log10", name = "Catch (clear) and SSB (solid)") + #, limits = c(0.5*min(plot_dat$value),NA)) +
  #   scale_fill_manual(name = "Species", values = params_uncalibrated@linecolour) +
  #   scale_alpha_manual(name = "Stat", values = c(0.5, 1), labels = c("Catch", "SSB")) +
  #   theme(
  #         legend.position = "none", legend.key = element_rect(fill = "white"))
    
 
ggplot(plot_dat)+
  geom_point(aes(x = w_inf, y = value, color = Species, shape = variable), size = 6, alpha = .8) +
  geom_point(data = plot_dat2, aes(x = w_inf, y = value*1562500, color = Species, shape = "averaged fishing mortality"), size = 6, alpha = .8)+
  geom_text_repel(data = filter(plot_dat,variable == "average SSB"), aes(x = w_inf, y = value, label = Species), hjust = 0, nudge_x = 0.05)+
  geom_line(aes(x = w_inf, y = value, color = Species)) +
  scale_y_continuous(name = "Catch and SSB", limits = c(0,NA),sec.axis = sec_axis(trans = ~./1562500)) +
  scale_x_continuous(name = "Asymptotic size (g)", trans = "log10") +
  scale_color_manual(name = "Species", values = params_uncalibrated@linecolour) +
  scale_shape_manual(name = "Data", values = c(16,4,17)) +
  theme(panel.background = element_blank(),
        panel.grid.minor = element_line(color = "gray"),
        legend.position = "bottom",legend.key = element_rect(fill = "white"))+
  guides(color = FALSE)
```

The first plot shows the fishing time-series per species while the second plot shows the average ssb and catch for the period 2014 - 2019.


RF: Dab and Gurnard have some issues, it comes from the data though, the calculations are correct

## Step 2. Guessing the maximum recruitment 

In this section you will:

- guess reasonable $R_{max}$ values which will stop out-competition from a few species

- take a look how it affects the density dependence of each species and their predicted yield



```{r step 2 - guessing coexistence, warning=F, message=F}
# let's start again and replace with the initial pre-calibration "guessed" Rmax 
params_guessed <- params_uncalibrated
# penalise the large species with higher density dependence
params_guessed@species_params$R_max <- params_guessed@resource_params$kappa*params_guessed@species_params$w_inf^-1
# and reduce erepro
# params_guessed@species_params$erepro <- 1e-3
params_guessed <- setParams(params_guessed)
# run with fishing
sim_guessed <- project(params_guessed, t_max = 100, effort =1)
plotSummary(sim_guessed, short = T)
```

The ecosystem looks way better. Saithe's largest individuals are having a hard time, but at least species coexist.

```{r step 2 - RDI/RDD}
  plot_dat <- as.data.frame(getRDD(sim_guessed@params)/getRDI(sim_guessed@params))
  plot_dat$species <- factor(rownames(plot_dat),sim_guessed@params@species_params$species)
  colnames(plot_dat)[1] <- "ratio"
  plot_dat$w_inf <- sim_guessed@params@species_params$w_inf
# ggplot(plot_dat) +
#   geom_bar(aes(x = species, y = ratio, fill = species),stat="identity") +
#   scale_fill_manual(name = "Species", values = sim_guessed@params@linecolour) +
#   scale_y_continuous(name = "density-independent / density-dependent reproduction rate", trans = "log10") +
#   scale_x_discrete(name = "Species") +
#    theme(panel.background = element_rect(fill = "white", color = "black"),
#         legend.position = "none")
ggplot(plot_dat)+
  geom_point(aes(x = w_inf, y = ratio, color = species), size = 6, alpha = .8) +
  geom_text(aes(x = w_inf, y = ratio, label = species), hjust = 0, nudge_x = 0.05)+
  # geom_line(aes(x = w_inf, y = value, color = Species)) +
  scale_y_continuous(name = "density-independent / density-dependent reproduction rate", limits = c(0,1)) +
  scale_x_continuous(name = "Asymptotic size (g)", trans = "log10") +
  scale_color_manual(name = "Species", values = params_uncalibrated@linecolour) +
  theme(panel.background = element_blank(),
                   panel.grid.minor = element_line(color = "gray"),
       legend.position = "none")
```


Is the physiological recruitment, $RDI$, much higher than the realised recruitment, $RDD$? High $RDI/RDD$ ratio indicates strong density dependence, meaning that the carrying capacity is controlling the population rather than predation or competition. Larger species often require more of this density dependent control than smaller ones. If $RDI/RDD$ is too high, the efficiency of reproduction ($erepro$) can be lowered to ensure species do not outcompete others or or over-resilient to fishing. The largest species that were the most limited by our new $R_{max}$ do not show a strong density dependence. The medium-sized species are the most affected here.




```{r step 2 - observed/predicted catches}
plotPredObsYield(sim_guessed,catchAvg$Catch_1419_tonnes) 
```


This plot show the ratio between observed yield (from the data loaded earlier) and the yield recreated in our model. Only one species matches at the moment (Sole). The size of the dot on the plot is proportional to the species' asymptotic size for clarity.


## Step 3. Calibrate the maximum recruitment

In this section you will:

- use a package that will calibrate $R_{max}$ per species

$Rmax$ will affect the relative biomass of each species (and, combined with the fishing parameters, the catches) by minimising the error between observed and estimated catches or biomasses. We could also include $\kappa$ in our estimation here (as in Blanchard et al 2104 & Spence et al 2016) but instead we will use the value that seemed OK in terms of feeding levels in the Rshiny app, roughly $log10(11.5)$. Same goes for $erepro$, a value of $1e-3$ seemed ok.

First let's set up a function running the model and outputing the difference between predicted catches (`getYield()`) and actual catches (`catchAvg`). `err` is the sum of squared errors between the two.

```{r step 3 - getError | function, include=F}
## the following getError function combines the steps of the optimisastion above - this time with the multispecies model and output the predicted size spectrum
## update below with project_steady and saving the state from each iteration
#RF the function takes a bunch of RMax and compare the theoretical catches versus data
getError <- function(vary,params,dat,env=state,data_type="catch", tol = 0.1,timetorun=10) {
  
  #env$params@species_params$R_max[]<-10^vary[1:12]
  params@species_params$R_max[]<-10^vary[1:12]
  
  params <- setParams(params)
  # run to steady state and update params
  # env$params<- projectToSteady(env$params, distance_func = distanceSSLogN,
  #                 tol = tol, t_max = 200,return_sim = F)
  params<- projectToSteady(params, distance_func = distanceSSLogN,
                   tol = tol, t_max = 200,return_sim = F)
 
  # create sim object 
   
  sim <- project(params, effort = 1, t_max = timetorun) #Change t_max to determine how many years the model runs for
  
  # 
  # sim <- project(env$params, effort = 1, t_max = timetorun) #Change t_max to determine how many years the model runs for
  # 
  # env$params <-sim@params
  # 
  
          ## what kind of data and output do we have?
          if (data_type=="SSB") {
          output <-getSSB(sim)[timetorun,]   #could change to getBiomass if using survey, also check units.
          }
  
          if (data_type=="catch") {
         output <-getYield(sim)[timetorun,]/1e6 
         #' using n . w . dw so g per year per volume (i.e. North Sea since kappa is set up this way). 
         #'The data are in tonnes per year so converting to tonnes.
          }
  
  pred <- log(output)
  dat  <- log(dat)
  # sum of squared errors, here on log-scale of predictions and data (could change this or use other error or likelihood options)
   discrep <- pred - dat
   discrep <- (sum(discrep^2))
  
  # can use a strong penalty on the error to ensure we reach a minimum of 10% of the data (biomass or catch) for each species
  # if(any(pred < 0.1*dat)) discrep <- discrep + 1e10
  
    return(discrep)
   }
```

```{r step 3 - getError | result}
# we need 12 Rmaxs, log10 scale
vary <- log10(params_guessed@species_params$R_max)
#vary<-runif(10,3,12) # or use completley made up values, same for each species test for effects of initial values
## set up the enviornment to keep the current state of the simulations 
state <- new.env(parent = emptyenv())
state$params <-  params_guessed
#catchAvg <-read.csv("data/time-averaged-catches.csv") # only that one is used at the moment | catches are estimated from fMatW
## test it
err<-getError(vary = vary, params = params_guessed, dat = catchAvg$Catch_1419_tonnes)
# err<-getError(vary,params,dat=rep(100,12),data_type="biomass")
err
```


Now, carry out the optimisation. There are several optimisation methods to choose from - we need to select the most robust one to share here. The R package optimParallel seems to be the most robust general R package and has replaced optim. Often this requires repeateing the proceure several times but the advantage of using parallel run is the speed compared to packages such as optimx.

This might take AWHILE. The output is saved as "optim_para_result" if you wish to skip this block.

```{r step 3 - optimisation, message = F, eval=F}
# change kappa and erepro based on shiny exploration, set up initial values based on "close to" equilibrium values from above sim
# params_steady already set to erepro = 0.001 and kappa = 10^11
params_optim <- params_guessed
vary <-  log10(params_optim@species_params$R_max)
# params_optim@resource_params$kappa<-3.2e11 # better kappa estimated from Rshiny
params_optim<-setParams(params_optim)
noCores <- detectCores() - 1 # keep a spare core
cl <- makeCluster(noCores, setup_timeout = 0.5)
setDefaultCluster(cl = cl)
clusterExport(cl, as.list(ls()))
clusterEvalQ(cl, {
  library(mizerExperimental)
  library(optimParallel)
})
tic()
optim_result <-optimParallel(par=vary,getError,params=params_optim, dat = catchAvg$Catch_1419_tonnes, method   ="L-BFGS-B",lower=c(rep(3,12)),upper= c(rep(15,12)),
                            parallel=list(loginfo=TRUE, forward=TRUE))
stopCluster(cl)
toc() # 80'' using 47 cores
saveRDS(optim_result,"optim_para_result.RDS")
```


```{r step3 - results, message=FALSE,warning=FALSE}
# if previous block was not evaluated
params_optim <- params_guessed
# params_optim@resource_params$kappa<-3.2e11 # old value (keeping it for now cause it works)
optim_result <- readRDS("optim_para_result.RDS")
# optim values:
params_optim@species_params$R_max <- 10^optim_result$par 
# set the param object 
params_optim <-setParams(params_optim)
sim_optim <- project(params_optim, effort = 1, t_max = 100, dt=0.1,initial_n = sim_guessed@n[100,,],initial_n_pp = sim_guessed@n_pp[100,])
saveRDS(sim_optim,"optim_para_sim.RDS")
plotSummary(sim_optim, short = T)
plotPredObsYield(sim_optim,catchAvg$Catch_1419_tonnes) 
```

Calibrating for best observed/predicted yield makes one species (Sprat) show signs of collapse. We need to look at other parameters to get the community to coexist again.



## Step 4. Calibrate the recruitment with eRepro.

In this section you will:

- Look at effect of $erepro$ on the reproductive outputs

- Check what impact erepro has on the $F_{msy}$



$eRepro$ represents the energy conversion efficiency between energy allocated to reproduction and the actual spawn. Lowering $erepro$ biologically means higher egg mortality rate or wasteful energy invested into gonads. For example, $eRepro$ is currently set to $0.01$ meaning, that for every one $g$ of mass allocated to reproduction, $0.1g$ will be used as spawn recruitment, the rest is lost. This coefficient is applied before any density-dependence requirement. Let's use Rshiny to see how varying $eRepro$ influences the ecosystem (to play around).


RF: shiny disabled for pdf


```{r step 4 r-shiny, eval = F, cache = F, include = F}
# Optional
# adjust Rmax and/or reproductive efficiency to examine whether improved steady state is achievable that way
library(shiny) # no need if runtime = shiny is in the YAML
# runApp("shiny-equilibrium")
# is there a way to save the final chosen values?
params_shiny <- params_optim
shinyApp(
  ui=fluidPage(
  
  # Application title
  titlePanel("North Sea Model Example"),
  
  fluidRow(
    column(4, wellPanel(
       sliderInput("kappa", "log10 Resource Carrying Capacity:", min = 8, max = 12, value = log10(params_optim@resource_params$kappa),
                   step = 0.1),
    #   sliderInput("Rmax", "log10 Maximum Recruitment:", min = 1, max = 12, value = 12,
    #              step = 0.1),
       sliderInput("erepro", "log10 Reproductive Efficiency:", min = -8, max = 1, value = log10(params_optim@species_params$erepro[1]),
                   step = 0.1)
          )),
    column(6,
           plotOutput("distPlot", width = 600, height = 600)
    ))
     
  
    
  ),
  
  server = function(input, output) {
   
  output$distPlot <- renderPlot({
    # set up params using values given, need check and change parameter values so units work in days units 
    params_shiny@species_params$erepro <- rep(10^input$erepro,12)
   # params@species_params$Rmax <- rep(10^input$Rmax,12)
    params_shiny <- setParams(params_shiny,kappa=10^input$kappa)
    # run without fishing
    sim_shiny <- project(params_shiny, effort = 1, t_max = 100)
    plot(sim_shiny)
     })
},
  options = list(height = 500)
)
```


Now let's take a look at the current $F_msy$ per species. To do so we need to set up an ecosystem with one gear per species and vary one gear intensity at a time while the other gear stay constant (default effort value of one).

```{r step 4 - Fmsy | data crunching, eval=F}
# Can be a bit long to not need to evaluate, just load the next block
plot_dat <- plotFmsy(params_optim, returnData = T, effortRes = 50)
saveRDS(plot_dat, "Fmsy.rds")
```


```{r step 4 - Fmsy | plots, warning=F}
plot_dat <- readRDS("Fmsy.rds")
ggplot(plot_dat) +
 geom_line(aes(x = effort , y = yield, color = species))+
  facet_wrap(species~., scales = "free") +
  scale_x_continuous(limits= c(0,1.5),name = "fishing mortality rate")+#, limits = c(1e10,NA))+
  scale_y_continuous(trans = "log10") +
  scale_color_manual(name = "Species", values = params_optim@linecolour) +
    theme(legend.position = "none", legend.key = element_rect(fill = "white"),
           panel.background = element_blank(), panel.grid.minor = element_line(color = "gray"),
          strip.background = element_blank())
```

Decreasing $eRepro$ is going to move the $MSY$ towards lower effort. Until now we had the same $eRepro$ for all species but we can input species-specific values to calibrate our recruitment. To do so, let's use another shiny app again. The next one allows you to change $eRepro$ one at a time and see the effect on the species' $F_{msy}$

RF: shiny disable for pdf


```{r step 4 - shiny Fmsy, eval = F, cache = F, include = F}
# Optional
# adjust Rmax and/or reproductive efficiency to examine whether improved steady state is achievable that way
library(shiny) # no need if runtime = shiny is in the YAML
# runApp("shiny-equilibrium")
# is there a way to save the final chosen values?
params_shiny <- params_optim
FmsyDat <- readRDS("Fmsy.rds")
shinyApp(
  ui=fluidPage(
  
  # Application title
  titlePanel("North Sea Fmsy"),
  
  fluidRow(
    column(4, wellPanel(
       # sliderInput("kappa", "log10 Resource Carrying Capacity:", min = 8, max = 12, value = log10(params_optim@resource_params$kappa),
       #             step = 0.1),
    #   sliderInput("Rmax", "log10 Maximum Recruitment:", min = 1, max = 12, value = 12,
    #              step = 0.1),
       sliderInput("erepro", "log10 Reproductive Efficiency:", min = -8, max = 1, value = log10(params_optim@species_params$erepro[1]),
                   step = 0.1),
sliderTextInput(
    inputId = "species",
    label = "Species name",
    choices = params_shiny@species_params$species,
    selected = params_shiny@species_params$species[1],
    grid = T
  ),
          )),
    
    column(6,
           plotOutput("distPlot", width = 600, height = 600)
    ))
     
  
    
  ),
  
  server = function(input, output) {
   
  output$distPlot <- renderPlot({
    # set up params using values given, need check and change parameter values so units work in days units 
    params_shiny@species_params$erepro[which(params_shiny@species_params$species == input$species)] <-10^(input$erepro)
    
    plotFmsy(params_shiny, speciesData = list(input$species,FmsyDat))
     })
},
  options = list(height = 500)
)
# add species slider to display jsut one species at a time
```


Let's see what our system looks like with species-specific $eRepro$

```{r step 4 - Fmsy with tweaked erepro, warning=F}
params_optim@species_params$erepro <- 10^(c(-1,-2,-4,-4,-3,-3,-2,-3,-2,-2,-2,-2))
params_optim2 <- setParams(params_optim)
sim_optim2 <- project(params_optim2, effort = 1, t_max = 100)
plotSummary(sim_optim2, short = T)
plotFmsy(params_optim2, effortRes = 30, returnData = F)
```

Some trade-offs here, relatively high $eRepro$ allows Sprat to coexist with the othe species but he $F_{msy}$ plot looks wrong. The same is true for Dab but I do not want to lower too much $eRepro$ either. Instead we are going to run the $R_{max}$ calibration again to see if we improved the yield match.

And now let's use this set of parameters to run the $R_{max}$ calibration again as in step 3.


```{r step 4 - rmax calibration crunch, message = F, eval=F}
# saving the last time step of sim_optim2 in the param object
params_optim2@initial_n <- sim_optim2@n[dim(sim_optim2@n)[1],,]
params_optim2@initial_n_pp <- sim_optim2@n_pp[dim(sim_optim2@n_pp)[1],]
params_calibration <- params_optim2
vary <-  log10(params_calibration@species_params$R_max)
params_calibration<-setParams(params_calibration)
noCores <- detectCores() - 1 # keep a spare core
cl <- makeCluster(noCores, setup_timeout = 0.5)
setDefaultCluster(cl = cl)
clusterExport(cl, as.list(ls()))
clusterEvalQ(cl, {
  library(mizerExperimental)
  library(optimParallel)
})
tic()
optim_result <-optimParallel(par=vary,getError,params=params_calibration, dat = catchAvg$Catch_1419_tonnes, 
                             method   ="L-BFGS-B",lower=c(rep(3,12)),upper= c(rep(15,12)),
                            parallel=list(loginfo=TRUE, forward=TRUE))
stopCluster(cl)
toc() # 80'' using 47 cores
saveRDS(optim_result,"optim_para_result2.RDS")
```

```{r step 4 - rmax calibration results, warning=F}
optim_result <- readRDS("optim_para_result2.RDS")
# if previous block not evaluated
params_optim2@initial_n <- sim_optim2@n[dim(sim_optim2@n)[1],,]
params_optim2@initial_n_pp <- sim_optim2@n_pp[dim(sim_optim2@n_pp)[1],]
params_calibration <- params_optim2
# optim values:
params_calibration@species_params$R_max <- 10^optim_result$par 
# set the param object 
params_calibration <-setParams(params_calibration)
sim_optim2 <- project(params_calibration, effort = 1, t_max = 100, dt=0.1, initial_n = params_calibration@initial_n ,
                      initial_n_pp = params_calibration@initial_n_pp)
saveRDS(sim_optim2,"optim_para_sim2.RDS")
plotSummary(sim_optim2, short = T)
plotPredObsYield(sim_optim2,catchAvg$Catch_1419_tonnes) 
```

Now the fit is less good overall but we only have 2 outliers, N.pout and Sandeel


## Step 5. Calibrating the growth

In this section you will:

- look at the growth curves of each species

- tweak the $\kappa$ parameter to adjust the growth curves


```{r step 5 - growth curves prior}
plotGrowthCurves2(sim_optim2, species_panel = T)
# Ken's comment: plotGrowthCurves. Legend should be “model”, “Observed”. Need to discuss about that since we need to change the plot in Mizer then
```

Most species have a growth curve similar to the expected von Bertalanffy growth curve, apart from Spart and Sandeel which have much slower growth than expected. $\kappa$, which is the carrying capacity of the background spectrum will affect the food availability and therefore the growth rate (more food means faster growth). However, changing $\kappa$ will affect all growth curves. If only a few species are off, we need to change the $\gamma$ parameter (search volume) per species.



```{r step 5 - feeding level}
plotFeedingLevel(sim_optim2, include_critical = F)
```

The feeding level here shows that the species reaching the highest size classes have a feeding level close to one, meaning that they feed to satiation. Let's look at their diet to check what makes them to full.

```{r step 5 - diets, warning=F}
# plotDiet2(sim_optim2, "Haddock")
# plotDiet2(sim_optim2, "Plaice")
# plotDiet2(sim_optim2, "Saithe")
plotDiet2(sim_optim2)
```

The diets look great as the fish feed on each other and do not rely too much on the background spectrum. We can still reduce the carrying capacity of the background spectrum even more, which could lower the feeding level of the species. Let's do this using the Rshiny app


RF: Rshiny section disabled for pdf

```{r step 5 r-shiny, eval = F, cache = F, include = F}
# Optional
# adjust Rmax and/or reproductive efficiency to examine whether improved steady state is achievable that way
library(shiny) # no need if runtime = shiny is in the YAML
# runApp("shiny-equilibrium")
# is there a way to save the final chosen values?
params_shiny <- sim_optim2@params
shinyApp(
  ui=fluidPage(
  
  # Application title
  titlePanel("North Sea Model Example"),
  
  fluidRow(
    column(4, wellPanel(
       sliderInput("kappa", "log10 Resource Carrying Capacity:", min = 8, max = 12, value = log10(params_shiny@resource_params$kappa),
                   step = 0.1),
    #   sliderInput("Rmax", "log10 Maximum Recruitment:", min = 1, max = 12, value = 12,
    #              step = 0.1),
       # sliderInput("erepro", "log10 Reproductive Efficiency:", min = -8, max = 1, value = -2,
       #             step = 0.1)
          )),
    column(6,
           plotOutput("distPlot", width = 600, height = 600)
    ))
     
  
    
  ),
  
  server = function(input, output) {
   
  output$distPlot <- renderPlot({
    # set up params using values given, need check and change parameter values so units work in days units 
    # params_shiny@species_params$erepro <- rep(10^input$erepro,12)
   # params@species_params$Rmax <- rep(10^input$Rmax,12)
    params_shiny <- setParams(params_shiny,kappa=10^input$kappa)
    # run without fishing
    sim_shiny <- project(params_shiny, effort = 1, t_max = 100)
    plot(sim_shiny)
     })
},
  options = list(height = 500)
)
```

Decreasing slightly $\kappa$ drives Sprat to extinction probably as it becomes the next choice of food when we reduce the size of the background spectrum (as we saw on the diet plots). So how can we reduce the feeding level of the largest predators while keeping Sprat alive?

```{r step 5 - kappa tweak, include = F}
# updating param parameters
params_optim2@resource_params$kappa<- 10^(11.5)
params_optim2 <-setParams(params_optim2)
sim_optim2 <- project(params_optim2, effort = 1, t_max = 100, dt=0.1)#,initial_n = sim_optim@n[100,,],initial_n_pp = sim_optim@n_pp[100,])
# saveRDS(sim_optim,"optim_para_sim2.RDS")
plotGrowthCurves(sim_optim2, species_panel = T)
# keep last time step of sim_optim saved to not have to get Sprat coming back from the deads every time
params_optim2@initial_n <- sim_optim2@n[dim(sim_optim2@n)[1],,]
params_optim2@initial_n_pp <- sim_optim2@n_pp[dim(sim_optim2@n_pp)[1],]
```




## Step 6. Calibrate food competition (rPP)

check what happens when we play with rPP

look at food level at larvea size and keep it constant by increasinf kappa when playing with rpp, the rest of the background will adjust itself from this one consaant size

look at resource mortality, choose favorite resource size, that size choose wich rpp give right plankton moratlity

## Step 7. Diagnostic plot



```{r step 7 - plots, message = F, warning=F, eval = F}
plotSummary(sim_optim)
plotPredObsYield(sim_optim,catchAvg$Catch_1419_tonnes) 
plotDiet2(sim_optim,"Cod") 
plotGrowthCurves(sim_optim, species_panel = T)
```

```{r step 7 - plotly, eval = F}
# interactive plots / won't be displayed in pdf
plotlyBiomass(sim_optim)
plotlySpectra(sim_optim)
plotlySpectra(sim_optim,power=2,total = T)
plotlyGrowthCurves(sim_optim,percentage = T) 
plotlyFeedingLevel(sim_optim) 
plotlyPredMort(sim_optim)
plotlyFMort(sim_optim)
# What would happen if we also parameterised the interaction matrix or beta and sigma?
```