
---
title: "StrathE2E2 User Manual (Version 3.0.0)"
author: |
  | Michael R. Heath (m.heath@strath.ac.uk) and Ian Thurlbeck (i.thurlbeck@strath.ac.uk)
  | University of Strathclyde, Department of Mathematics and Statistics,
  | Glasgow, Scotland, UK
date: "`r format(Sys.Date(), '%B %d, %Y')`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{StrathE2E2 User Manual}
  %\VignetteEngine{knitr::rmarkdown}
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(out.extra = '', fig.pos = 'H', collapse = TRUE, comment = NA)
```

\


# Summary

* **StrathE2E2** is an **R** package for modelling the 'big-picture', whole ecosystem effects of hydrodynamics, temperature, nutrient additions, and fishing 
on continental shelf marine food webs.

* **StrathE2E2** has two linked parts - a fishing fleet model and an ecology model. The fishing model integrates harvesting, discarding and seabed 
disturbance rates across a range of gears and passes the results into the ecology model.

* The ecology model is a network of coupled ordinary differential equations representing the rates of change in organic detritus, dissolved inorganic 
nutrient, and coarse guilds of living biomass spanning microbes to megafauna. The equations include representations of feeding, metabolism, reproduction, active migrations, advection and mixing. Environmental driving data include temperature, irradiance, hydrodynamics, and nutrient inputs from rivers, atmosphere and ocean boundaries.

* The package includes functions for parameter optimization, global sensitivity analysis, Monte Carlo estimation of credible intervals for model outputs.

* A fully developed and documented implementation of StrathE2E for the North Sea in included with the package.

\


# Documentation

**To see an overview of the package, launch an R-session and then:**    
```{r launchoverview1, eval = FALSE}
library(StrathE2E2)     # Load the package
help(StrathE2E2)        # Launch the overview document in a browser window
```

\

[**Cheatsheet** - quick-reference guide to StrathE2E2 functions](StrathE2E2_Cheatsheet.pdf)  

\

[**Technical Manual** - documentation on the structure of input and output R-objects and files](StrathE2E2_Technical_Manual.pdf)  

\

## See also the following supporting documents:
[Origin of the model (previous versions)](StrathE2E2_Model_origin.html)  
[Overview of the model concepts and design](StrathE2E2_Model_overview.html)  
[Ecology model description](StrathE2E2_Ecology_model.pdf)  
[Fishing fleet model description](StrathE2E2_Fishing_fleet_model.pdf)  
[Parameter optimization, sensitivity and Monte Carlo methodology](StrathE2E2_Optimization_sensitivity_Monte_Carlo_methodology.pdf)  
[StrathE2E implementation for the North Sea](StrathE2E2_North_Sea_model.pdf)

\


## Funding
Development of the StrathE2E2 package was supported by the [UK NERC Marine Ecosystems Research Programme (NE/L003120/1)](https://www.marine-ecosystems.org.uk/Home)    
Other projects which have supported aspects of the model development are:  
[EU-Horizon 2020 (DISCARDLESS project)](http://www.discardless.eu/)   
[Fisheries Innovation Scotland (project FIS003)](https://fiscot.org/wp-content/uploads/2019/06/FIS003.pdf).

\


## User community
We plan to develop a repository of open access implementations of StrathE2E for shelf sea regions. Users developing implementations 
which they are willing to share are encouraged to contact us to discuss inclusion in the repository.

\


## How to cite the package:
```{r launchoverview2, eval = TRUE}
citation("StrathE2E2")
```

\

\


# Contents

## Quick start

## 1 Package structure

## 2 Getting started and model house-keeping
### 2.1 Installing the package
### 2.2 Folder structure for model input and output files
### 2.3 Listing and copying model folders
#### 2.3.1 Listing the available models
#### 2.3.2 Copying a named model
### 2.4 Downloading parameter documentation

## 3 Basic model operations
### 3.1 Reading a model configuration
### 3.2 Running a model
### 3.3 Setting initial conditions 
### 3.4 Extracting harvest ratios from a model run
### 3.5 Plotting time-series of model outputs

## 4 More advanced use of the model
### 4.1 Writing your own code to create model scenarios
### 4.2 Writing your own code to work with the model output

## 5 Parameter estimation
### 5.1 Background
### 5.2 Implementation of the Metropolis-Hastings algorithm
### 5.3 Target data for optimization
### 5.4 Estimating of ecology model parameters
### 5.5 Estimating of harvest ratio scaling parameters
### 5.6 Estimating of fishing gear activity rates
### 5.7 Diagnostics for optimization processes
### 5.8 Calculating initial values for harvest ratio scaling parameters

## 6 Sensitivity and Monte Carlo analysis
### 6.1 Background
#### 6.1.1 Sensitivity analysis
#### 6.1.2 Monte Carlo analysis
### 6.2 Performing a sensitivity analysis
### 6.3 Performing a Monte Carlo analysis
### 6.4 Parallelisation of sensitivity and Monte Carlo analyses
#### 6.4.1 Merging data generated on parallel processors
#### 6.4.2 Processing data generated on parallel processors
### 6.5 Diagnostics for sensitivity and Monte Carlo analyses

## 7 Comparing model runs and observations
### 7.1 Comparing model runs with observations
### 7.2 Comparing baseline and scenario model runs
#### 7.2.1 Comparing baseline and scenario model results using boxplots
#### 7.2.2 Comparing baseline and scenario model results using tornado barplots

## 8 Generating fisheries yield curves
### 8.1 Generating yield data
### 8.2 Plotting yield data

## 9 Plotting functions for visualizing model inputs
### 9.1 Visualizing environmetal drivers
### 9.2 Visualizing fishing-related drivers

## 10 Plotting functions for visualizing model outputs
### 10.1 Plotting annual cycles of ecology model variables
### 10.2 Plotting annual cycles of active migration fluxes
### 10.3 Plotting distribitions of catches by gear and guild
### 10.4 Plotting mean trophic level and omnivory indices
### 10.5 Plotting zonal distributions of annual average biomass densities

## References

\

\


# Quick start
```{r quickstart, fig.width=7, fig.height=6, fig.align='center'}
library(StrathE2E2)				# Load the package
model <- e2e_read("North_Sea", "1970-1999")	# Read the 1970-1999 version of the North Sea model
results <- e2e_run(model,nyears=3)		# Run the model for 3 years
e2e_plot_ts(model,results,"ECO")	        # Plot time series of the masses of the state variables
```

\

# 1 Package structure

The package contains 8 families of functions:

* Model house-keeping (**Table 1**)
* Basic model operations (**Table 2**)
* Parameter estimation (**Table 3**)
* Sensitivity and Monte Carlo analyses (**Table 4**)
* Comparing model runs and observations (**Table 5**)
* Fishery yield curves (**Table 6**)
* Visualizing model inputs (**Table 7**)
* Visualizing model outputs (**Table 8**)

\

**Table 1**. Model house-keeping

| Function name     | Description                                                                         |
| ----------------- | ----------------------------------------------------------------------------------- |
| e2e_ls()          | List the available models                                                           |
| e2e_copy()        | Copy models between workspaces                                                      |
| e2e_get_parmdoc() | Download parameter documentation as a dataframe                                     |

\

\

**Table 2**. Basic model operations

| Function name       | Description                                                                  |
| ------------------- | ---------------------------------------------------------------------------- |
| e2e_read()          | Read a model setup                                                           |
| e2e_run()           | Run StrathE2E for a prescribed number of years with a given setup            |
| e2e_extract_start() | Create a new initial values file from the end-state of a model run           |
| e2e_extract_hr()    | Extract the values of harvest ratios generated by the fleet model            |
| e2e_plot_ts()       | Plot time series of model outputs for the full duration of a model run       |

\

\

**Table 3**. Parameter estimation
 
| Function name              | Description                                                                         |
| -------------------------- | ----------------------------------------------------------------------------------- |
| e2e_optimize_eco()         | Optimize ecology model parameters                                                   |
| e2e_optimize_hr()          | Optimize harvest ratio scaling parameters                                           |
| e2e_optimize_act()         | Optimize fishing gear activity rates                                                |
| e2e_plot_opt_diagnostics() | Plot diagnostic data for optimization runs                                          |
| e2e_calculate_hrscale()    | Calculate initial values of parameters linking effort to harvest ratios             |


\

\

**Table 4**. Senstivity and Monte Carlo analyses

| Function name           | Description                                                           |
| ----------------------- | --------------------------------------------------------------------- |
| e2e_run_sens()          | Run a global parameters sensitivity analysis with a given setup       |
| e2e_run_mc()            | Run a Monte Carlo analysis with a given setup                         |
| e2e_merge_sens_mc()     | Merge parallel processing files from sensitivity or Monte Carlo runs  |
| e2e_process_sens_mc()   | Post-process raw output data from sensitivity or Monte Carlo runs     |
| e2e_plot_sens_mc()      | Plot diagnostic results from sensitivity or Monte Carlo runs          |

\

\
 
**Table 5**. Compare model runs and observations

| Function name              | Description                                                        |
| -------------------------- | ------------------------------------------------------------------ |
| e2e_compare_obs()          | Box-plot comparisons between observations and model outputs        |
| e2e_compare_runs_box()     | Box-plot comparisons between two different model runs              |
| e2e_compare_runs_bar()     | Tornado bar-plots comparing two different model runs               |

\

\
 
**Table 6**. Fishery yield analysis

| Function name              | Description                                                        |
| -------------------------- | ------------------------------------------------------------------ |
| e2e_run_ycurve()           | Perform a set of model runs to generate fishery yield curve data   |
| e2e_plot_ycurve()          | Plot fishery yield curve data                                      |

\

\

**Table 7**. Visualize model inputs

| Function name              | Description                                                        |
| -------------------------- | ------------------------------------------------------------------ |
| e2e_plot_edrivers()        | Plot a climatological year of environmental driving data           |
| e2e_plot_fdrivers()        | Plot distributions of fishery-related driving data                 |

\

\

**Table 8**. Visualize model outputs from the final year of a run

| Function name              | Description                                                              |
| -------------------------- | ------------------------------------------------------------------------ |
| e2e_plot_eco()             | Plot annual cycles of ecology model variables                            |
| e2e_plot_catch()           | Plot annual landings and discards                                        |
| e2e_plot_migration()       | Plot annual cycles of active migration fluxes                            |
| e2e_plot_trophic()         | Plot mean trophic level and omnivory indices                             |
| e2e_plot_biomass()         | Plot zonal distributions of annual average biomass densities             |

\

\

# 2 Getting started and model house-keeping

## 2.1 Installing the package

Download the installation zip file for your operating system and install into your R version. Note that **Strathe2E2** requires installation 
of the following additional packages which are available from any CRAN site:

* deSolve

* NetIndices

\

To start a session, begin by loading the package:
```{r loadpackage}
library(StrathE2E2)
```

\


## 2.2 Folder structure for model input and output files

The data which define a model (parameters and driving data) are held in a specific folder structure for each implementation. The location of the 
North Sea demonstration model data within the folder structure of the R-package is shown in **Figure 1**.

\

![](images/folder_structure_figure_1.PNG)

**Figure 1**. Locations of the definition data folders for the North Sea demonstration model and its two variants in the structure of the 
R-package.

\

The folder structure within the model package is read-only, so users developing their own models will need to create an editable folder structure in 
their own accessible workspace to hold their model configuration data. In addition, users of the North Sea demonstration model wishing to use the more 
advanced optimization functions in the package will need to create a copy of the demonstration models in their own workspace, because these functions 
require to write data back into the **/Parameters** sub-folders on completion of the optimization. This can be done using using any file manager or the 
**e2e_copy()** function provided with the package. The structure expected for any user-defined model folders is shown in **Figure 2**.
 
\

![](images/folder_structure_figure_2.PNG)

**Figure 2**. Organisation of model definition data folders in any user-defined work-space. Italicised folder names are free for user specification.

\

The model outputs R-list objects which are held in memory and available for user-analysis. In addition, many of the contents of the list objects are output 
to comma-separated-variable (csv) files in a folder structure in the current user workspace (**Figure 3**). If a results folder does not exists the package will auto-create one:

\

![](images/folder_structure_figure_3.PNG)

**Figure 3**. Organisation of model results folders.

\

\

## 2.3 Listing and copying model folders

### 2.3.1 Listing the available models **e2e_ls()**

List the available models in a designated workspace.

**Table 9**. Arguments of the **e2e_ls()** function

| Argument       | Description                                                                               |
| -------------- | ----------------------------------------------------------------------------------------- |
| path           | Path to users model folder, otherwise read the list of models supplied with the package   |

\

```{r listmodels}
library("StrathE2E2")
e2e_ls()
```

\

\

### 2.3.2 Copying a named model

Make a copy of a named model/variant and documentation in a user defined workspace.

**Table 10**. Arguments of the **e2e_ls()** function

| Argument       | Description                                                                               |
| -------------- | ----------------------------------------------------------------------------------------- |
| model.name     | Name of a model to copy  |
| model.variant  | Name of a model variant to copy
| dest.path      | Destination folder to write a copy of the model, which will be created if necessary. The default is "Models" in the current folder |
| user.path      | Source folder o fthemodel to be copied. If the source is one ofthe North Sea models provided with the package then user.path=""

\

```{r copymodels, eval = FALSE}
library("StrathE2E2")
e2e_copy(" Modelname " , " Modelvariant ", dest.path= " mymodels ", user.path= " source " )
#
# The following call would copy one of the example models in he package to a user location:
e2e_copy("North_Sea" , "1970-1999", dest.path= " mymodels " )
```

\

Copies of the North Sea model make a useful template for development of new models.
 
\

\

## 2.4 Downloading parameter documentation **e2e_get_parmdoc()**

Get documentation on the parameters inthe model in a dataframe format

**Table 11**. Arguments of the **e2e_ls()** function

| Argument     | Description                                                             |
| ------------ | ----------------------------------------------------------------------- |
| id           | Integer value denoting the class of parameter details to be downloaded  |

\

The full range of inputs to the model includes both numeric values which are constant over the during a run, and values which vary over time during a run
according to an externally provided schedule. We refer to the latter as 'drivers', but they are all parameters nontheless.
The constant parameters are divided into categories according to whether they are associated with the ecology model, the fishing fleet model, and whether 
they are available for optimization or not.

\

This function downloads documentation on the full range of parameters in the form of a dataframe which can be used to form labels in user-generated code. 
The parameters are identified by a numeric value linked to 6 different classes of parameter. 
Choose from: 0 = fitted ecology, 1 = fixed ecology, 2 = fishing fleet, 3 = harvest ratios, 4 = environmental drivers, 5 = physical configuration, 9 = all (default = 9).

\

Note that the function does not supply the numeric values of parameters for any given model setup. These are available in the model definition object generated by the function e2e_read().

\

```{r getparmdoc}
parm_list <- e2e_get_parmdoc(0)       # Get documentation on the fitted ecology parameters
parm_list <- e2e_get_parmdoc(2)       # Get documentation on the fishing fleet model parameters
parm_list <- e2e_get_parmdoc()        # Get documentation all parameters
head(parm_list)
```

\

\

# 3 Basic model operations

## 3.1 Reading a model configuration **e2e_read()**

Read all the configuration, driving and parameter files for a designaded model and compile into a list object.

**Table 12**. Arguments of the **e2e_read()** function

| Argument       | Description                                                                                                                                             |
| -------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------- |
| model.name     | Name of model to read (no default)                                                                                                                      |
| model.variant  | Variant of the model to be read (no default)                                                                                                            |
| model.ident    | Text appended to output csv file names (e.g. OFFSHORE_model_annualresults-TEXT.csv instead of just OFFSHORE_model_annualresults.csv) (default="base") |
| model.subdir   | Sub directory of the main results folder for storing results (default = a folder which will be auto-created using the model name)                       |
| user.path      | Path to users top level model folder from where the model setup will be loaded (default = load from the package folder extdata/Models)                  |
| quiet          | Set to FALSE to see on-screen information on where model is loaded from and what files are read (default=TRUE)                                          |

\

The **e2e_read()** function is the foundation for all StrathE2E2 model runs. This function loads all the data required to configure a model into an R-list object. 
The function is pointed to each of the required input files by the file **MODEL_SETUP.csv** which resides in the **/ModelName/VariantName** folder for the relevant 
model. **MODEL_SETUP.csv** points to specific files in the **/Parameters**, **/Target_data** and **/Driving_data** subfolders.

\

The default **user.path** for **e2e_read()** is to load one of the variants of the North Sea model supplied with the package. The North Sea models are 
configured so that the initial conditions which are loaded for a default run launch the system in a stationary state - that is, the model produces 
repeating annual cycles of output given repeating annual cycles of the time-dependent driving data (e.g. temperature, hydrodynamics etc).

\

```{r loadmodel}
library("StrathE2E2")
model <- e2e_read("North_Sea", "1970-1999")
```

The object **model** is a list containing all of the input data and parameters gathered from the csv files in the folder structure:

```{r viewmodelobject}
str(model,max.level=2)                    # View the top 2 levels of input list object
```

\

Full documentation of the input csv files and the contents of the list object produced by **e2e_read()** is provided in the [Technical Manual](StrathE2E2_Technical_Manual.pdf).

\


## 3.2 Running a model **e2e_run()**

Run the StrathE2E model with a given configuration.

**Table 13**. Arguments of the **e2e_run()** function

| Argument       | Description                                                                            |
| -------------- | -------------------------------------------------------------------------------------- |
| model          | R-list object defining the model configuration compiled by the e2e_read() function   |
| nyears         | Number of years (integer) to run the model (default=20)                                |
| quiet          | Set to TRUE to quieten status outputs during the model run (default=TRUE)              |
| csv.output     | Set to FALSE to disable writing of csv output files (default=TRUE)                     |

\

The function **e2e_run()** runs a single instance of themodel, and generates another list object containing all of the outputs. In addition a range 
of derived outputs for the final year of the run are exported as csv files to a **/results** folder in the users local workspace.

\


```{r runmodel}
results <- e2e_run(model, nyears=3)         # Run the model for 3 years
str(results,max.level=1)                    # View the high-level contents of the output list object
```

\

Full documentation of the output csv files and the contents of the list object produced by **e2e_run()** is provided in the [Technical Manual](StrathE2E2_Technical_Manual.pdf).
\


## 3.3 Setting initial conditions **e2e_extract_start()**

Extract the values of all state variables at the end of a model run and format for use as new initial conditions.

**Table 14**. Arguments of the **e2e_extract_start()** function

| Argument   | Description                                                                         |
| ---------- | ----------------------------------------------------------------------------------- |
| model      | R-list model object for the model run generated by the function **e2e_read()**      |
| results    | R-list results object generated by the function **e2e_run()**                       |
| csv.output | Logical. If FALSE then disable writing of csv file (default = TRUE)                 |

\

The function **e2e_extract_start()** saves the state of the model at the end of a run (using the **e2e_run()** function) for use as initial conditions in 
future runs. This enables, for example, the model to be run for a long time to attain a stationary state, and then restarted in that state.

\

Initial conditions for a model run are held in the **/Parameters** folder of the Model/Variant path specified in the **e2e_read()** 
function call used to define a run. By default, the function attampts to write the model end-state file back to this **/ModelName/VariantName/Parameters** folder.
However, the package folders are read-only so if **e2e_read()** has been specified to load an internally provided **Model/Variant** then 
the output will revert to the currently specified results folder instead. To fix this, copy the required package model to a user
workspace using the **e2e_copy()** function and re-run.

\

The internal North Sea models are both supplied with initial conditions which are consistent with a stationary state, so they should both show repeating annual cycles of
output from time = 0, unless the parameters or driving conditions are altered to create some new environmental scenario.

\

```{r extractendcon, eval=FALSE}
model<-e2e_read("North_Sea", "2003-2013")                             # Read the 2003-2013 internal model
results <- e2e_run(model, nyears=3)                                   # Run the model for 3 years
new_start_data <- e2e_extract_start(model,results,csv.output=FALSE)   # Extract the model endstate as a dataframe only
# This will produce the new data only as a dataframe
head(new_start_data)
#
#
# Copy the 2003-2013 version of the North Sea model supplied with the package into a user workspace (Windows OS):
e2e_copy("North_Sea", "2003-2013",
            dest.path="C:/Users/username/Documents/Models")
# Load the 2003-2013 version of the North Sea model from the user workspace:
model<-e2e_read(model.name="North_Sea",
                  model.variant="2003-2013",
                  model.ident="TEST",
                  user.path="C:/Users/username/Documents/Models")
results <- e2e_run(model, nyears=3)                    # Run the model for 3 years
new_start_data <- e2e_extract_start(model,results)     # Extract the model endstate as new initial conditions file
# This will write to the user model /Parameters folder)
```

\

\

## 3.4 Extracting harvest ratios from a model run **e2e_extract_hr()**

Extract the values of harvest ratios generated by the fleet model.

**Table 15**. Arguments of the **e2e_extract_hr()** function

| Argument   | Description                                                                         |
| ---------- | ----------------------------------------------------------------------------------- |
| model      | R-list model object for the model run generated by the function **e2e_read()**      |
| results    | R-list results object generated by the function **e2e_run()**                       |
| csv.output | Logical. If FALSE then disable writing of csv file (default = TRUE)                 |

\

Harvest ratios for each guild in each zone (proportion of biomass captured per day) are a key input to the ecology model. They are generated by the fleet model 
based on inputs of fishing gear activity rates, distribution in space, and catching power for each guild.

\

The function **e2e_extract_hr()** returns a dataframe, and optionally a csv file, of harvest ratios in each zone and the domain as a whole, for each guild in the ecology model. These are extracted 
from the results object generated by running the **e2e_run()** function.

\

```{r extractharvestratios, eval=TRUE}
# Load the 2003-2013 version of the North Sea model from the user workspace:
model<-e2e_read(model.name="North_Sea",
                  model.variant="2003-2013")
results <- e2e_run(model, nyears=3)                                   # Run the model for 3 years
harvest_ratio_data <- e2e_extract_hr(model,results,csv.output=FALSE)  # Extract the harvest ratio data
harvest_ratio_data
```

\

\

## 3.5 Ploting time-series of model outputs **e2e_plot_ts()**

Plot ecology model state variables or fishery catches for the full duration of a model run.

**Table 16**. Arguments of the function for plotting full-length time-series data from model outputs

| Argument  | Description                                                                           |
| --------- | ------------------------------------------------------------------------------------- |
| model     | R-list model object for the model run generated by the function **e2e_read()**      |
| results   | R-list results object generated by the function **e2e_run()**                       |
| selection | Text string from a list identifying whether ecology state variable or fishery catches are to be plotted. Select from: "ECO", "CATCH", default = "ECO". Remember to include the phrase within "" quotes. |

\

The function **e2e_plot_ts()** plots a multi-panel page of time series plots of either a) daily ecology state variable values aggregated over the inshore and offshore zones of the domain and over sediment classes, 
or b) annual landings and discards by zone for each model guild, for the full duration of a model run. Currently the masses of macrophytes, corpses and discards 
are not included in the state varibale plots due to space constraints.

\

Be warned that if the run is more than about 10 years then the plot becomes extremely compressed and messy. It is not intended to
be of publication quality. The intention is to provide a quick-look diagnostic of trends in the state variables. This is useful to assess
whether the model is close to stationary state or not.

\

Units of the plotted state varaiables are mMN i.e. mass in the model domain without any scaling to zone-area or layer thickness. Similarly, units of catches 
are mMN.y^-1^ from the whole model domain without any scaling to zone-area. The assumed area of the whole model domain is 1 m^2^.

\

Selection of either ecology state variable or catches to plot is by a function argument.

\

A separate set of functions is provided for plotting more detailed visualizations of just the final year from a model run, e.g. **e2e_plot_eco()**, **e2e_plot_catch()**.

\

```{r plot_ts_01, fig.width=7, fig.height=6, fig.align='center'}
model<-e2e_read(model.name="North_Sea",
                  model.variant="1970-1999")
results <- e2e_run(model, nyears=3)
e2e_plot_ts(model, results,selection="ECO")
```

**Figure 4**. Time-series plot of ecology model output, units mMN in the model domain as a whole (1 m^-2^ sea surface).

\

```{r plot_ts_02, fig.width=7, fig.height=6, fig.align='center'}
e2e_plot_ts(model, results,selection="CATCH")
```

**Figure 5**. Time-series plot of catch data output, units mMN.y^-1^ in the model domain as a whole (1m^-2^ sea surface).

\

All of the plotting functions provide with the package produce output to a screen window. To divert output to a graphic file wrap the plotting function call 
inside graphics device open and close statememts as follows:

```{r plot_ts_02_tofile, eval = FALSE}
pdf("plot.pdf",width=7,height=6)                    # or jpeg("plot.jpg"), png("plot.PNG")
e2e_plot_ts(model, results,selection="CATCH")
dev.off()
```

\

\

# 4 More advanced use of the model

## 4.1 Writing your own code to create model scenarios

Typical use of the model might involve

* comparing a baseline run with a scenario run involving some changes in driving data (e.g. different temperature conditions or different activity rates of selected gears),

* conducting an analysis of the sensitivity to systematic changes in driving data (e.g. incrememts in river nutrient inputs over some range of values).

\

Users have two options for configuring such experiments

* manually edit the .csv input files as required and rerun the model,

* use coding to alter the required elements of the R list object created by the **e2e_read()** function, assigning a unique identifier to the outputs from each run. This is the most efficient way to create model scenarios

\

The scope for developing fishing and/or environmemtal scenarios by coding to modify a default model object loaded by **e2e_read()** is limitless, but requires 
and understanding of the structure and content of the list object produced by the **e2e_read()** function. Use the **str()** function in R to explore the list and 
its sub-components, and consult the [Technical Manual](StrathE2E2_Technical_Manual.pdf).

\

```{r viewmodelobject2}
model<-e2e_read("North_Sea", "2003-2013")
str(model,max.level=2)                    # View the top 2 levels of input list object
```

\

To mine deeper into the model object, for example, to show the fleet model inputs:

```{r viewfleetmodelinput}
model<-e2e_read("North_Sea", "2003-2013")
model$data$fleet.model
```
\

Below are two examples of code to configure and run scenario cases of the the 2003-2013 North Sea model provided with the package. 
The first example runs 2003-2013 North Sea model as a baseline, and then adds 2 deg-C to the temperatures in all spatial zones and 
reruns the model. The second example iterates through seven levels of demersal fish harvest ratio ranging from 0 to 3-times the 
value in the 1970-1999 model (in steps of 0.5-times). Note that the changes in harvest ratio are not accompanied by any changes in 
fishing gear activity. The package includes a dedicated function for performing this task and generating plots (**e2e_run_ycurve()**.

\

NOTE - The model is run for only 10 years in Example 1 to create the output shown here - in an actual working example it would be 
recommended to run for at least 40 years in each case to ensure that the model is close to a stationary state by the end of each run.

\

Model scenario Example 1

```{r runmodelscenario1-1, fig.width=7, fig.height=6, fig.align='center'}
# Example of code to run a baseline case of the North Sea model with 2003-2013 conditions, and then edit the
# model list object to create a scenario run with the temperature in all zones increased by 2 deg-C
#--------------------------------------------------------
# Read the embedded North Sea 2003-2013 model and assign an identifier for the results
baseTemp_model<-e2e_read("North_Sea", "2003-2013",model.ident="baseTemp")
# Run the model for 10 years and save the results to a named list object
baseTemp_results<-e2e_run(baseTemp_model,nyears=10)
# Visualize the output from the run (should show a repeating annual cycle with no trend)
e2e_plot_ts(baseTemp_model,baseTemp_results,"ECO")
```

**Figure 6**. Baseline model run for Example 1 - 2003-2013 conditions in the North Sea

\

```{r runmodelscenario1-2, fig.width=7, fig.height=6, fig.align='center'}
T_rise<-2    # temperature increase to add to all temperatures in the model drivers
# copy the baseline model list object to a new model list
baseTemp_plusTC_model<-baseTemp_model
# add temperature increase to upper layer offshore temperatures
baseTemp_plusTC_model$data$physics.drivers$so_temp <- baseTemp_model$data$physics.drivers$so_temp+T_rise
# add temperature increase to inshore temperatures
baseTemp_plusTC_model$data$physics.drivers$si_temp <- baseTemp_model$data$physics.drivers$si_temp+T_rise
# add temperature increase to lower layer offshore temperatures
baseTemp_plusTC_model$data$physics.drivers$d_temp  <- baseTemp_model$data$physics.drivers$d_temp+T_rise
# Assign a unique identifier for the .csv outputs
baseTemp_plusTC_model$setup$model.ident <- "baseTemp_plusTC"
# Run the model for 10 years and save the results to a named list object
baseTemp_plusTC_results<-e2e_run(baseTemp_plusTC_model,nyears=10)
# Visualize the output from the run (should show trends in outputs due to change in T)
e2e_plot_ts(baseTemp_plusTC_model,baseTemp_plusTC_results,"ECO")
```

**Figure 7**. Scenario model run for Example 1 - 2003-2013 plus 2^o^C conditions in the North Sea.

\

Model scenario Example 2

```{r runmodelscenario2, eval = FALSE}
# Example of code to loop through a set of seven levels of demersal fish harvest ratio ranging from 0 to 3-times
# the baseline value for the 1970-1999 North Sea model. The .csv outputs for each level are saved to unique
# filenames but the results list-object is discarded. The baseline case is the third level (HScale=1). Each level
# is run for 40 years. 
#--------------------------------------------------------
model<-e2e_read("North_Sea", "1970-1999")  		# Read the embedded North Sea 1970-1999 model
for(i in 1:7) {                              		# Start loop through seven model scenarios
  HScale <- (i-1)*0.5                          		# HScale becomes 0 to 3 in steps of 0.5
      # Assign a unique identifier for each set of .csv outputs
  model$setup$model.ident <- paste("Dem_HR_scale_",HScale,sep="")
      # Set the demersal fish harvest ratio multipler to HScale
  model$data$fleet.model$HRscale_vector_multiplier[2] <- HScale
  results<-e2e_run(model,nyears=40)			 # Run the model, save the .csv results, overwrite any previous list-object
  cat("Harvest Ratio scaling = ",HScale,"\n")            # Print a screen message to monitor progress
  e2e_plot_ts(model,results,"ECO")	  	         # Visualize the output from each run
}							 # End loop through scenarios
```

\

\

## 4.2 Writing your own code to work with the model output

The first step to working with the model outputs is understanding their structure. This is described in the following section, but more details are available in a 
[Technical Manual](StrathE2E2_Technical_Manual.pdf).

\

The outputs from running the basic model are an R list object (**results** in the examples above) and a set of csv files. The list object 
contains 7 primary data structures. The csv file outputs are generated in the **/results** folder, either in the user home workspace by default, or in a workspace 
defined by the arguments of the **e2e_read()** function. Each csv file has a generic file name plus an alpha-numeric identifier (**model.ident**) assigned by the user 
as an argument of the **e2e_read()** function (default = "base"). In the following descriptions this is represented by "TEXT". The csv files replicate elements of the 
data contained in the R list object generated by **e2e_run()**. Their purpose
is to preserve the most useful outputs for future use e.g. in comparisons of baseline and scenario model runs.

\

An abbreviated view of the contents of the results list object is shown below. 


```{r runmodel3}
model<-e2e_read("North_Sea", "2003-2013")
results <- e2e_run(model, nyears=3)         # Run the model for 3 years
str(results,max.level=1)                    # View the high-level contents of the output list object
```

\

\

A brief description of the contents of each of the elements of the results object, and correspondence with csv output files is shown below. Full details are provided in the [Technical Manual](StrathE2E2_Technical_Manual.pdf).

\

The object **build** is a list containing a copy of the model input data imported from the **model** object created by **e2e_read()**. This is included in 
the **results** object since some of these data are required for post-processing of results, e.g. to convert state variable mass to area or volume densities (see examples later):

```{r resultsobjectdata1}
str(results$build,max.level=1) 
```

\

The object **output** is a dataframe containing the raw daily-interval output from the ODE solving computations on the masses of all of the state varaiables. 
This dataframe is the raw material for all subsequent post-processing of the model results. The example below views just the top and bottom few rows, and the 
first 6 columns of the dataframe. 
The prefix x_ in front of the detritus and other terms signifies a sediment property, with the extension _s1, _s2, or _s3 indicating inshore sediment habitat, and 
_d1, _d2, or _d3 indicating offshore sediment habitat. For water column components the extension _so, _d, and _si indicate offshore upper and lower, and inshore zones/layers.

```{r resultsobjectdata2}
head(results$output[,1:6])
tail(results$output[,1:6])
```

\

The object **aggregates** is a list containing vectors of 124 derived model variable which are aggregations of selected columns in the raw model output in the
**output** object. Aggregations are created by summing, for example, across inshore and offshore values of the same food web component to produce a whole-domain value.
Each numeric vector has the same length as the number of rows in dataframe **results\$output**, in other words a value for each daily output step.
In addition, the **aggregates** list contains six single numeric values of zonal layer thicknesses and sediment porosities needed to convert aggregated masses into
mass densities (m^-2^ or m^-3^).

\

The example below views the names of just the first few variables in the **aggregates** list, extracts two variables and makes plot of their time-series. More details are provided in the [Technical Manual](StrathE2E2_Technical_Manual.pdf).

```{r resultsobjectdata3}
# List the names of the first 40 vectors in the results$aggregates list:
names(results$aggregates[1:40])
```

```{r plotagg, fig.width=6, fig.height=4, fig.align='center'}
# Plot the time series of the sedient nitrate mass in the inshore and offshore zones
# i.e. aggregated across the sediment types in each zone.
# Note that this is the mass not the porewater concentration so uncorrected for
# zonal area or sediment porosity or layer thickness.
plot(seq(1,nrow(results$output)), results$aggregates$x_nitrate_o, type="l",col="black",
     ylim=( c( 0,(1.5*max(results$aggregates$x_nitrate_o)) ) ),
     xlab="Days", ylab="Sediment nitrate mass (mMN)")
lines(seq(1,nrow(results$output)), results$aggregates$x_nitrate_i, type="l",col="red")
legend("topleft",  bg="transparent", c("offshore","inshore"), lty=c(1,1),col=c("black","red"))
```

**Figure 8**. Plot of the time series of sediment nitrate mass aggregated over the offshore and inshore zones.

\

The object **fleet.output** is a list of dataframes containing outputs from the fleet model which are either piped into the ecology model, or saved internally for use 
in post-processing the model outputrs - in particular for disaggregating the total catch of each good web guild into the fractions attributable to each fishing gear. 
This post processing is carried out automatically by the **e2e_run()** function (see below) so under most circumatances the user would not require access to these 
internal data:

```{r resultsobjectdata4}
str(results$fleet.output,max.level=1)
```

\

The object **total.annual.catch** is a list comprising two dataframes - **offshore_annual_group_gear_land_disc** and 
**inshore_annual_group_gear_land_disc**. These dataframes contain the annual landings and discards of each harvestable guild in the model 
for each year of the model run (units mMN.y^-1^, from a model domain of 1m^2^ sea surface area). These data correspond to the files 
**model_offshore_annual_landings_discards-TEXT.csv** and **model_inshore_annual_landings_discards-TEXT.csv** which are generated in the **/results** folder. 

```{r resultsobjectdata5}
str(results$total.annual.catch,max.level=1)
```

\

The object **annual.catch.by.gear** is a list comprising two dataframes - **offshore_annual_group_land_disc** and 
**inshore_annual_group_land_disc**. These dataframes contain the same data as in **total.annual.catch** but disaggregataed into the 
fractions attributable to the (up to) 12 gears. The data are not output to csv files.

```{r resultsobjectdata6}
str(results$annual.catch.by.gear,max.level=1)
```

\

The object **final.year.output** is a list of arrays and dataframes containing a wide range of properties derived from the final year of the model run. The majority of 
these are also output as csv files in the **/results** folder.

```{r resultsobjectdata7}
str(results$final.year.output,max.level=1)
```

\

The first six entries in the **final.year.output** list object are arrays which all have the same format - the catch, landings or discard mass of each guild (rows) by each gear (columns) in a given spatial zone, over the 
final year (mMN.y^-1^ from a model domain of 1m^2^ surface area). These data correspond with csv files in the **/results** folder, for example: 
**final.year.output\$inshore_catchmat** becomes **INSHORE_catchcomposition_by_gear-TEXT.csv**. For example, the first two columns of the array containing inshore catch data are given by:

```{r resultsobjectdata8}
# List the first six entries
str(results$final.year.output[1:6],max.level=1)
#
# View the first two column of the inshore_catchmat array
results$final.year.output$inshore_catchmat[,1:2]
```

\

The object **final.year.output\$monthly_averages** is a dataframe of monthly averaged nutrient and plankton data (units: mMN.m^-3^, except for chlorophyll: mg.m^-3^) corresponding to the observational data available in 
the **/Target_data** folder of model stup data. Chlorophyll is estimated assuming Redfield C:N (molar ratio 106:16), and a carbon:chlorophyll weight ratio of 20. The data correspnd to the csv output file 
**model_monthlyresults-TEXT.csv**. For example, the first two columns of the array are given by:

```{r resultsobjectdata9}
results$final.year.output$monthly.averages[,1:2]
```

\

The next four entries in the **final.year.output** list object are dataframes containing annual averages, maxima and minima of each state variable in the inshore zone (mMN), 
and a set of annually integratred flux terms (mMN.y^-1^) during the final year. The fluxes include both physical flows due to advection, mixing, and inflows and outflows, and ecological fluxes between state variables. 
The food web fluxes are not calculated for the inshore zone, but are provided in other data frames for the model domain as a whole.  In addition the dataframe contains the zonal area proportions and layer thicknesses 
needed to convert the data to mass or flux densities (m^-2^ or m^-3^). The formats of the average, maxima and minima files are the same. The data correspond to csv files output to the */results* folder with names: 
**INSHORE_model_anav_biomass-TEXT.csv**, **INSHORE_model_maximum_biomass-TEXT.csv**, **INSHORE_model_minimum_biomass-TEXT.csv**, **INSHORE_model_annualresults-TEXT.csv**. 
The dataframes contain 3 columns: the numeric values, units and a descriptive text field.
For example, list these dataframes in the list, and view the top few rows of the annual average dataframe, and the annual fluxes dataframe:

```{r resultsobjectdata10}
str(results$final.year.output[8:11],max.level=1)
#
head(results$final.year.output$mass_results_inshore)
#
head(results$final.year.output$annual_flux_results_inshore)
```

\

The next four entries in the **final.year.output** list object are dataframes containing annual averages, maxima and minima of each state variable, and annual fluxes in 
the offshore zone (mMN) during the final year, exactly corresponding to the previous four entries for the inshore zone. As for the inshore zone, the data are output to 
corresponding csv files.

```{r resultsobjectdata11}
str(results$final.year.output[12:15],max.level=1)
```

\

The next four entries in the **final.year.output** list object are dataframes containing annual averages, maxima and minima of each state variable, and annual fluxes for 
the model domain as a whole (mMN) during the final year, exactly corresponding to the previous four entries for the offshore zone, except that in this case all of the 
flows between state variables in the food web are evaluated. In following outputs these flows are assembled into the flow matrix format required by the NetIndices package 
to calculate network indices. As for the offshore zone, the data are output to corresponding csv files.

```{r resultsobjectdata12}
str(results$final.year.output[16:19],max.level=1)
```

\

The next three entries in the **final.year.output** list object are dataframes containing the flow matrices required for input to the **NetIndices** package, and the outputs from NetIndices. 
Two version of the flow matrix are included - one which includes spawning and recruitment flows - from adults to larvae and larvae to adults, and one in which these flows are excluded. 
In each case the rows names of the matrix indicate the "flow-from", and the columns indicatethe "flow-to". 
All three of these data sets are exported to csv files in the **/results** folder: **flow_matrix_all_fluxes-TEXT.csv**, **flow_matrix_excl_spawn_recruit-TEXT.csv**, and **Network_indices_output-TEXT.csv**. 

```{r resultsobjectdata13}
str(results$final.year.output[20:22],max.level=1)
#
#Show the first 5 rows and columns of the matrix of all flows:
results$final.year.output$flow_matrix_all_fluxes[1:5,1:5]
#
#Show the first few rows of the dataframe containing outputs from the NetIndices package:
head(results$final.year.output$NetworkIndexResults)
```

\

The next two entries in the **final.year.output** list object are dataframes containing annual observational data on the state of the ecosystem, and the monthly observational data. These are 
exact copies of the input data from the **/Target_data** folder of input files, and are included in the **results** object merely for convenience so that they are easily available
for comparison with model outputs:

```{r resultsobjectdata14}
str(results$final.year.output[23:24],max.level=1)
```

\

The final three elements of the **final.year.output** list object are associated with the calulation of the likelihood of the observed annual target data on the state of the ecosystem, 
given the model drivers and parameters. The first is a single numeric value - the overall likelihood of the observed target data. The second is a dataframe of the partial likelihoods 
for each aspect of the observational dataset. Finally, a dataframe of the derived model outouts for the final year of the run, exactly corresponding to the set of observational measures.
The two dataframes are also output as csv files to the **/results** folder with names **model_likelihood_results-TEXT.csv** and **model_target_annualresults_plus_chi-TEXT.csv**.

```{r resultsobjectdata15}
str(results$final.year.output[25:27],max.level=1)
#
#Show the overall likelihood value:
results$final.year.output$annual_obj
#
#Show the first few rows of the partial likelihood data:
head(results$final.year.output$partial_chi)
#
#Show the first few rows of the dataframe containing the model outputs corresponding to the annual observational target data:
head(results$final.year.output$opt_results)
```

\

While the **results** object generated by the **e2e_run()** function contains a wide range of derived model outputs, users may also need to generate their own. 
The following examples illustrate how to extract and process the raw output held in the **results\$output** dataframe.

\

For a set of results created using the **e2e_read()** and **e2e_run()** functions, the raw data from the model run are contained in the object 
dataframe **results\$output**. The first row of the **output** dataframe corresponds to time=0, i.e. the initial conditions passed to the differential equation solver 
at the start of the run. A calendar year of output corresponds to a time interval of 360 days, i.e. 361 rows of output data (time = 0:360), so the total number of 
rows of output will be ((360 * nyears)+1), where nyears is the specified length of the model run in years.

Assuming that the **output** dataframe has been extracted from the results object into a new dataframe **out** as follows:

```{r processresults1}	
out <- results$output
```

Data for any one year specified by Y (Y <= nyears) can be extracted by the R statement:

```{r processresults2}	
# All the output corresponding to year Y
Y <- 3                    # select year 3
yearY_data <- out[(((Y-1)*360)+1) : ((Y*360)+1),]
#
#List the first few rows and columns of these extracted data
head(yearY_data[,1:6])
```

\

```{r processresults3}	
# All the output for a given state variable (denoted by "column_name") in year Y:
Y <- 3                    # select year 3
column_name<- "herb_o"    # Select offshore omnivorous zooplankton
col2get <- which(colnames(out)==column_name)
column_data_yearY <- out[(((Y-1)*360)+1) : ((Y*360)+1),col2get]
```

\

```{r processresults4}	
# Calculate the annual average in year Y of the variable extracted in the previous step:
	mean(column_data_yearY)
```

\

The flux terms in the **output** dataframe are cumulative over the duration of the model run, beginning at 0 at time = 0. Hence, the 
time-series of instantaneous flux rates (d^-1^) is the sequence of increments between successive rows of output. In R this can be generated by, for example:

```{r processresults5, fig.width=6, fig.height=5, fig.align='center'}	
column_name<- "phytgrossprod_o"    		   # Select for example offshore gross production by phytoplankton
col2get <- which(colnames(out)==column_name)
J <- nrow(out)					   # number of rows of data in the "out" dataframe
temp <- rep(0,J)				   # creates a vector to hold temporary data
temp[1: (J-1)] <- out[2 : J, col2get]		   # copies data from the "out" into "temp" offset by 1 time increment
temp[J] <- 2*(out[J,col2get]) - out[J-1,col2get]   # fills in the terminal value of the "temp" vector by extrapolation
rate <- temp - out[,col2get]			   # "rate" is a vector of the increments between successive time steps
#
#Plot the original cumulative data and the derived rate
par(mfrow=c(2,1))
par(mar=c(4,4,1,1))
plot(seq(1,length(out[,col2get])),out[,col2get],type="l",xlab="Days",ylab="Cumulative flux")
plot(seq(1,length(rate)),rate,type="l",xlab="Days",ylab="Daily flux")
```

**Figure 9**. Upper panel: time series of the cumulative phytoplankton gross production as output from the model. Lower panel: same data converted into daily gross production rates.

\

To obtain the annual integral of a flux variable over a given year, the R statement would be:

```{r processresults6}
Y <- 3                                                  # select year 3
column_name<- "phytgrossprod_o"    		        # Select for example offshore gross production by phytoplankton
col2get <- which(colnames(out)==column_name)
out[((Y*360)+1),col2get] - out[(((Y-1)*360)+1),col2get] # integrated flux of the variable denoted by "column_name" over year Y
```

\

All of the terms in the **output** dataframe are in mass units (mMN or mMC) in the model domain which is is scaled to a sea-surface area of 1 m^2^. 
To express the outputs as area densities in the inshore or offshore zones (e.g. mMN.m^-2^) or layer concentrations (e.g. mMN.m^-3^), the mass values need 
to be re-scaled to the area-proportion of the relevant zone and/or layer thickness. In addition, terms relating to sediment or sediment porewater 
require to be scaled by the area-proportions of seabed habitats, sediment layer thicknesses and sediment porosity. These scaling parameters are 
accessible to the user as elements of the object **build\$model.parameters** which forms part of the list-object generated by a model run. Examples of 
R statements for converting mass outputs into area-densities and layer concentrations are shown below. 

\

```{r processresults7}
# Extract the relevant area and volumetric parameters from the $build$model.parameters object
inshore_area 	     <- as.numeric(results$build$model.parameters["shallowprop"])
offshore_upper_thick <- as.numeric(results$build$model.parameters["thik_so"])
offshore_lower_thick <- as.numeric(results$build$model.parameters["thik_d"])
inshore_sed1_area    <- as.numeric(results$build$model.parameters["area_s1"])
inshore_sed1_thick   <- as.numeric(results$build$model.parameters["thik_x_s1"])
inshore_sed1_poros   <- as.numeric(results$build$model.parameters["porosity_s1"])
#
# Time series of suspended detritus and bacteria concentration in the offshore upper layer (mMN.m-3)
detr_bact_conc_so<- out$detritus_so/((1-inshore_area)*offshore_upper_thick)
#
# Time series of the area-density of omnivorous zooplankton in the inshore zone (mMN.m-2)
oz_dens_i<- out$herb_i/inshore_area		
#
# Time series of the depth averaged concentration of omnivorous zooplankton in the offshore zone (mMN.m-3)
oz_conc_o<- out$herb_o /((1-inshore_area)*(offshore_upper_thick+offshore_lower_thick))
#
# Time series of the porewater concentration of ammonia in inshore sediment habitat 1 (mMN.m-3):
xamm_conc_s1<- out$x_ammonia_s1/(inshore_sed1_area*inshore_sed1_thick*inshore_sed1_poros)
#
# Time series of labile detritus & bacterial in inshore sediment habitat 1 as a % of sediment dry weight (%gN.g-1):
xdet_conc_s1<- 100*(out$xR_detritus_s1*(14/1000))/(inshore_sed1_area*inshore_sed1_thick*(1-inshore_sed1_poros)*2650000 )
# (nitrogen atomic weight = 14 g.mole-1; dry sediment density = quartz density = 2.65 x 1e6 g.m-3)	
```

\

\

# 5 Parameter estimation

## 5.1 Background

Three groups of parameters in the combined ecology and fishing fleet model system might be considered as candidates for estimation, depending 
on the availability of external data:

* Ecology model parameters
* Scaling parameters linking fishing effort to harvest ratios
* Fishing gear activity rates

\

It is clearly unrealistic to expect to be able to optimize all of these sets of parameters concurrently since there is a risk that some may 
be confounded. So, the R package provides several options for parameter estimation, each of which optimizes one of these groups 
of parameters, given knowledge on the other two.

\

The parameter estimation schemes are configured as a stationary state fitting method, though there is no intrinsic 
reason why they could not be adapted for time-series fitting. For the most part, the methodology involves likelihood estimation and a simulated annealing scheme 
using the Metropolis-Hastings iterative algorithm (Bertsimas & Tsitsiklis, 1993; Cerny, 1985; Kirkpatrick et al., 1983).

\


## 5.2 Implementation of the Metropolis-Hastings algorithm

The Metropolis-Hastings algorithm implemented with StrathE2E2 automates the acceptance or rejection of proposed new randomly generated sets of 
model parameters (**Figure 10**). The objective is to maximize the likelihood of the set of observed (target) indices given a vector of model parameter 
values. Starting with an initial informed-guess parameter vector, the procedure iterates through model runs, and at each iteration the 
parameters are independently 'jiggled' by proposing new values drawn from gaussian distributions of fixed coefficient of variation (standard deviation 
divided by the mean) around the current accepted values. There are some special-case parameters in StrathE2E2 which are bounded by known limits, and the preference parameters which require to be 
renormalized to sum to 1.0 for each predator after the addition of random 'jiggles'.

\

![](images/metropolis_hastings_flow_diagram.PNG)

**Figure 10**. Flow diagram illustrating the operation of the simulated annealing scheme for parameter optimization.

\

After running the model to a stationary state with the new parameter vector, and extracting simulated versions of the ($i$) observed ecosystem state 
indices ($1 < i < I$) (where the observed indices each have a standard deviation $\sigma_i$), an overall error term ($\chi_{\theta^*}^{2}$) is calculated 
for the model as follows: 

$$\chi_{\theta^*}^{2} = \frac{\sum_{i=1}^{i=I} \frac{(observed_{i} - simulated_{\theta^{*},i} )^{2}}{2\sigma_{i}^{2} }}{I} $$

\

The likelihood of the overall set of indices $P(observation|\theta^{*})$ is then $exp(-\chi_{\theta^*}^{2}$).

\

The change in performance of the model due to the jiggling of the proposed parameter set is measured by ($\chi_{\theta^*}^{2} - \chi_{\theta_{current}}^{2}$).
A simple 'hill-climbing' scheme would accept the proposed parameter set $\theta^*$ as a new version of the parameter set $\theta_{current}$ for the 
next iteration only if the likelihood ratio $LR = exp(\chi_{\theta^*}^{2} - \chi_{\theta_{current}}^{2}) > 1$.
Over many iterations the accepted parameter set should then migrate towards a maximum likelihood fit. However, such 
schemes are notorious for becoming trapped in local optima. Simulated annealing attempts to mitigate this risk by accepting a proportion of 
instances of proposed sets which produce worse results than current accepted (i.e.  $LR < 1$), thereby exploring a wider range of the parameter space. 
The process is an analogy of the way in which crystals form in a liquid metal as it cools. Taking $T[k]$ to represent the 'temperature' of the 
system at iteration $k$, the probability of accepting a proposal as the new current is equal to the modified likelihood ratio:
$$LR_m = exp((\chi_{\theta^*}^{2} - \chi_{\theta_{current}}^{2})/T[k])$$

\

The value of $T[k]$ is allowed to 'cool' as the iterations progress (increasing $k$). As $T[k]$ decreases towards 0 the system assumes a simple 
hill-climbing mode (i.e. only vectors $\theta^*$ producing likelihood ratios > 1 are accepted).  In the early stages, when $T(k) > 0$, a proportion 
of $\theta^*$ producing likelihood ratios < 1 are accepted. Practically, after each model run $\theta[k]$ = $\theta^*$ only if $LR_m \ge \lambda$, 
where $\lambda$ is drawn from a uniform random distribution between 0  and 1, otherwise $\theta[k]$ = $\theta[k-1]$ for the next iteration. 
At each iteration, both the currently accepted and the proposed parameter vectors are saved, along with the resulting likelihood of the 
observations given the parameters.

\

The sequence of (P(observation|$\theta^*$)[k]) should converge to a maximum value representing the best attainable fit between the model and the observations 
given the model structure and driving data. However, the rate at which convergence is achieved, and potentially the converged vector itself, 
are dependent on the cooling schedule and the magnitude of the jiggling coefficient of variation (cv) used to generate $\theta[k]$. A geometrically 
decreasing temperature is applied ($T[k] = r.T[k-1]$), and after experimentation with the value of r and cv, the  finally adopted values for 
fitting StrathE2E2 models to achieve consistency of convergence within a sensible number of iterations (<10000) are: r = 0.975, cv = 0.005. 
The best-fit value of $\theta[k]$ is regarded as being that remaining when no new parameter vectors are accepted within 200 iterations.

\


## 5.3 Target data for optimization

The target data for optimizing a particular model should all be located in the **/Models/Variant/Target_data** folder of the model definition. The individual files names
must all contain a standard text string plus a free text string for identification of different versions. The particular version required to be loaded for a given
optimization run should be identified inthe **MODEL_SETUP.csv** file in the **/ModelName/VariantName** folder.

\

In most instances, the target data to which the model is optimized during the simulated annealing process are
a set of up to 84 indices of the real-world state of the ecosystem and 
their associated uncertainties (standard deviations). This file is **named annual_target_data\*.csv** (where \* is a free text field for identification).
The ecosystem state indices include annual average biomasses or concentrations, annual integrated production rates of 
food web components, annual integrated dietary intakes, fishery landings and discards, denitrification rate, and other properties. The file also contains
text fields including a description, units, geographic region and time period of applicability, and the origins of the data that were used to generate the index.
For each index, the code calculates an equivalent value for the final year of each model run as the basis for the likelihood estimation. Individual indices can 
be included or excluded from the likelihood calculation by means of a switch (value 1 or 0) in the input file.

\

One of the availble methods for estimating fishing gear activity rates where these are unknown relies on target data concerning the harvest ratios on 
each guild in each zone of the model. These are contained in a file **fishing_fleet_target_harvest_ratios\*.csv** (where \* is a free text field for identification).

\

Finally, calculation of effort to harvest ratio scaling parameter values when activity, fishing power and harvest ratios are all known, requires the target data on harvest ratios. 
to be located in the file **known_regional_harvest_ratios\*.csv** (where \* is a free text field for identification).

\


## 5.4 Estimating of ecology model parameters **e2e_optimize_eco()**

Optimize StrathE2E ecology model parameters to maximize the likelihood of observed ecosystem target data.

**Table 17**. Arguments of the function **e2e_optimize_eco()**.

| Argument          | Description                                                                                                                                                                                                                               |
| ----------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| model             | R-list object generated by the e2e_read() function which defines the model configuration                                                                                                                                                  |
| nyears            | Number of years to run the model in each iteration (default=40)                                                                                                                                                                           |
| n_iter            | Number of iterations of the model (default=500)                                                                                                                                                                                           |
| start_temperature | Initial value of the simulated annealing temperature parameter (default=0.5). Suggested values in the range 0.0005 - 5. Higher values increase the probability of rejecting parameter combinations producing an improvement in likelihood |
| cooling           | Rate at which the simulated annealing temperature declines with iterations (default=0.975). Suggested values in the range 0.9 - 0.985                                                                                                     |
| toppredlock       | Logical. If TRUE then locks-down the uptake parameters of the birds pinnipeds and cetaceans as these are hard to fit alongside the other parameters (default=TRUE)                                                                        |
| quiet             | Logical. If TRUE then suppress informational messages at the start of each iteration (default=TRUE)                                                                                                                                       |
| csv.output        |Logical. If FALSE then disable writing of csv output files - useful for testing (default=TRUE)                                                                                                                                             | 
| runtime.plot      |Logical. If FALSE then disable runtime plotting of the progress of the run - useful for testing (default=TRUE)                                                                                                                             |


The function **e2e_optimize_eco()** launches a StrathE2E simulated annealing process to find the set of ecology model parameters 
producing the maximum likelihood of observed target data on the state of the ecosystem, 
given specified environmental driving data and fishing fleet parameters. 

\

Certain parameters of the ecology model are excluded from the optimization process because they are either confounded with other parameters 
in some way or there is a clear route to establishing values from external evidence (**Table 18**). All other parameters, defining the prey 
preferences of all the living components of the model, the uptake and mortality rate coefficients, and the rate coefficients for the 
biogeochemical processes, are eligible for optimization. A subset of the feeding rate parameters for the top-predators can be optionally 
excluded from the optimization process if reasonable values can be independently estimated. Separate coefficients of variation can be set 
for different classes of parameters (preference parameters, maximum uptake rates, half saturation coefficients, microbial rates, density 
dependent mortality rates) via a setup text file.

\

**Table 18**. Background to the parameters of the ecological model which require to be fixed from external evidence and are excluded from the simulated annealing process. 
These parameters are contained in two files **fixed_parameters_for_consumer_groups\*.csv**  and **fixed_parameters_miscellaneous\*.csv** (where * is a free text identifier)
in the **/Parameters** folder of each model definition.

| Parameter type                                                                            | Description                                                                                                                                                 |
| ----------------------------------------------------------------------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------- |
| Assimilation efficiencies for each living guild in the model                              | Proportion of the ingested mass of food that contributes to new body tissue after subtracting defecation and the metabolic costs of digestion and synthesis |
| Biomass loss rates due to temperature-dependent metabolism for each living resource guild | Proportion of biomass lost to ammonia per day due to non-feeding related metabolism at a given Q10 reference temperature                                    |
| Q10 values for temperature dependent processes and the Q10 reference temperature          | Separate Q10 values for autotrophic uptake of nutrient, heterotrophic feeding, and heterotrophic metabolism based on literature data                        |
| Light intensity required to saturate autotrophic nutrient uptake                          | Light saturation intensity for nutrient uptake, estimated from survey of field and laboratory experiments                                                   |
| Annual weight specific fecundities of fish and benthos guilds                             | Guild-level values of the annual proportion of biomass which is shed as eggs                                                                                |
| Spawning start and end dates for fish and benthos guilds                                  | Day-number (assuming 360-day year) of the onset and end of spawning to be based on such empirical data as are available                                     |
| Recruitment start and end dates for fish and benthos guilds                               | Day-number (assuming a 360-day year) of the onset and end of recruitment from the larval stage to the main guild biomass                                    |
| Extra-domain biomass of migratory fish and the proportion invading the domain each year   | Biomass of the wider stock contributing to the transient population in the model domain, extracted from fisheries stock assessments                         |
| Start and end dates for the annual invasion and emigration of migratory fish              | Day-number (assuming a 360-day year) of the onset and end of immigration and emigration phases                                                              |
| Harvestable biomass density threshold for each resource guild                             | A limit for each guild below which the harvestable species are assumed to be exhausted. Values based on data from surveys, fisheries and stock assessments  |
| Minimum inedible biomass of carnivorous zooplankton                                       | A minimum edible threshold is set to ensure that the guild as a whole cannot be extirpated by predation                                                     |

\

Parameters elligible for fitting are contained in the files **fitted_parameters_preference_matrix\*.csv**, **fitted_parameters_uptake_and_mortality_rates\*.csv**, and 
**fitted_parameters_microbiology_and_others\*.csv** (where \* is a free text identifier) in the **/Parameters** folder of each model definition. At the end of an optimization 
run the code exports a new set of these three files containing the accepted values for the maximum likelihood model, back to the **/Parameters** folder with a new text
identifier which is set by the **model.ident** argument in the **e2e_read()** function.

\

Configuration of a simulated annealing run to estimate ecology model parameters must be preceded by a **e2e_read()** function call to load the model setup
which includes the initial values for all of the parameters to be fitted. It is important to note that the **user.path** argument in the **e2e_read()** function 
call needs to point to a user workspace folder, not the default
North Sea model provided with the package. This is because the annealing function needs write-access to the model **/Parameters** folder,
but the **/extdata/Models** folder in the package installation is read-only. To use the annealing function on the North Sea model, use the **e2e_copy()** 
function to make a copy of the North Sea model in the user workspace.

\

The coefficients of variation for 'jiggling' the parameter values during an optimization run are held in the file **annealing_SD_ecology\*.csv** in the 
**/Parameters/Parameter_SD_control** folder of each model definition. These can be changed in real-time during a run, to enable exploration of the effects of parameter 
search patterns.

\

The function produces a real-time graphical summary of the progress of the fitting procedure, displaying the likelihoods of the proposed and accepted 
parameter sets at each iteration. The Y-axis (likelihood) range of the real time plot can be varied during the run by editing the setup file **annealing_SD_ecology.csv** at any time during the run.

\

At the end of the procedure the function returns a list object containing the histories of proposed and accepted parameters and the final accepted parameter 
values. Optionally (by default), new versions of the three ecology model **fitted_parameter\*.csv** files are exported
to the folder **/Parameters** of the model version, with a user defined identifier specified by the **model.ident** argument
in the **e2e_read()** function. The histories of proposed and accepted parameter combinations
are saved as csv files in the results folder.

\

In order to use the new fitted parameter values in a subsequent run of the StrathE2E model (using the **e2e_run()** function) 
it will be necessary to edit the **MODEL_SETUP.csv** file in the relevant **/Models/variant** folder to point to the new files.

\

```{r annealingcode1, eval = FALSE}
# Copy the 2003-2013 version of the North Sea model supplied with the package into a user workspace (Windows OS):
	e2e_copy("North_Sea", "2003-2013",
	dest.path="C:/Users/username/Documents/Models")
#
# Load the 2003-2013 version of the North Sea model from the user workspace:
	model<-e2e_read(model.name="North_Sea",
                          model.variant="2003-2013",
                          model.ident="TEST",
                          user.path="C:/Users/username/Documents/Models")
#
# Quick demonstration of the annealing function in operation:
e2e_optimize_eco(model, nyears=5, n_iter=10, start_temperature=0.5)
#
# More realistic configuration would be as follows :
e2e_optimize_eco(model, nyears=50, n_iter=1000, start_temperature=1)
# (WARNING - this will take about 26 hours to run)
#
```

\

\

## 5.5 Estimating of harvest ratio scaling parameters **e2e_optimize_hr()**

Optimize StrathE2E harvest ratio multipliers to maximize the likelihood of observed ecosystem target data.

**Table 19**. Arguments of the function **find_harvest_ratio_mult()**.

| Argument          | Description                                                                                                                                                                                                                               |
| ----------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| model             | R-list object generated by the e2e_read() function which defines the model configuration                                                                                                                                                  |
| nyears            | Number of years to run the model in each iteration (default=40)                                                                                                                                                                           |
| n_iter            | Number of iterations of the model (default=500)                                                                                                                                                                                           |
| start_temperature | Initial value of the simulated annealing temperature parameter (default=0.5). Suggested values in the range 0.0005 - 5. Higher values increase the probability of rejecting parameter combinations producing an improvement in likelihood |
| cooling           | Rate at which the simulated annealing temperature declines with iterations (default=0.975). Suggested values in the range 0.9 - 0.985                                                                                                     |
| quiet             | Logical. If TRUE then suppress informational messages at the start of each iteration (default=TRUE)                                                                                                                                       |
| csv.output        |Logical. If FALSE then disable writing of csv output files - useful for testing (default=TRUE)                                                                                                                                             | 
| runtime.plot      |Logical. If FALSE then disable runtime plotting of the progress of the run - useful for testing (default=TRUE)                                                                                                                             |


\

The function **e2e_optimize_hr()** launches a StrathE2E simulated annealing process to find the set of fishing fleet model harvest ratio multipliers
producing the maximum likelihood of observed target data on the state of the ecosystem, given specified environmental driving data and ecology model parameters.

\

Proposed values of the scaling parameters are contained in the configuration file for the fishing fleet model. In addition, the model setup 
includes a separate csv file (**harvest_ratio_multiplier.csv**) containing multipliers to be applied to each of these scaling parameters. 
These are useful when conducting scenario experiments, but would ordinarily be set to 1.0. The simulated annealing code for optimizing the 
scaling values iterates the multiplier values to seek the combination producing the best fit of the overall model to the set of target data. 
At the end of the process, the code saves the history of proposed and accepted multiplier values to the current **/results** folder, and the final accepted 
version of the multipliers back to the **/ModelName/VariantName/Parameters** folder with the model.itent text identifier assigned in the **e2e_read()** function call which must precede
the run. If deemed worthy, these multipliers can then be manually applied to the corresponding parameter values in the fishing fleet configuration to produce an updated 
version. The effort-harvest ratio multipliers associated with birds, pinnipeds and cetaceans are likely to be poorly constrained by observational data
and may need to be set independently. Nevertheless thay are included in the process here.

\

Configuration of a simulated annealing run to estimate harvest ratio multipliers must be preceded by a **e2e_read()** function call to load the model setup
which includes the initial values of the parameters to be fitted. It is important to note that the **user.path** argument in the **e2e_read()** function 
call needs to point to a user workspace folder, not the default
North Sea model provided with the package. This is because the annealing function needs write-access to the model **/Parameters** folder,
but the **/extdata/Models** folder in the package installation is read-only. To use the annealing function on the North Sea model, use the **e2e_copy()** 
function to make a copy of the North Sea model in the user workspace.

\

The coefficient of variation for 'jiggling' the parameter values during an optimization run is held in the file **annealing_SD_fishing\*.csv** in the 
**/Parameters/Parameter_SD_control** folder of each model definition. This can be changed in real-time during a run, to enable exploration of the effects of parameter 
search patterns.

\

The function produces a real-time graphical summary of the progress of the fitting procedure, displaying the likelihoods of the proposed and accepted 
parameter sets at each iteration. The Y-axis (likelihood) range of the real time plot can be varied during the run by editing the setup file **annealing_SD_fishing.csv** at any time during the run.

\

At the end of the procedure the function returns a list object containing the histories of proposed and accepted parameters and the final accepted parameter 
values. Optionally (by default), new versions of the harvest ratio multiplier files is exported
to the folder **/Parameters** of the model version, with a user defined identifier specified by the **model.ident** argument
in the **e2e_read()** function. The histories of proposed and accepted parameter combinations
are saved as csv files in the results folder.

\

In order to use the new fitted parameter values in a subsequent run of the StrathE2E model (using the **e2e_run()** function) 
it will be necessary to edit the **MODEL_SETUP.csv** file in the relevant **/Models/variant** folder to point to the new files.

\

Configuration of a simulated annealing run to estimate the effort to harvest ratio scaling parameters by optimizing to the ecosystem state target data must be 
preceded by a **e2e_read()** function call to load the model setup. This will include the initial values for the harvest ratio multipliers which should ideally (though 
not necessarily be all set to 1.0 from the start.

\

It is important to note that the user.path argument in the **e2e_read()** function call needs to point to a user workspace folder, not the default
North Sea model provided with the package. This is because the annealing function needs write-access to the model **/Parameters** folder,
but the **/extdata/Models** folder in the package installation is read-only. To use the annealing function on the North Sea model, use the **e2e_copy()** 
function to make a copy of the North Sea model in the user workspace.

\

The coefficients of variation for jiggling the harvest ratio parameters can be varied in real-time during the run by editing the file **annealing_SD_fishing.csv** 
in the folder **/Parameters/Parameter_SD_control** of the model version.

\

The function produces a real-time graphical summary of the progress of the fitting procedure, displaying the likelihoods of the proposed and accepted 
parameter sets at each iteration. The Y-axis (likelihood) range of the real time plot can be varied during the run by editing the setup file **annealing_SD_fishing.csv** at any time during the run.

\

At the end of the procedure new versions of the harvest ratio multiplier file is exported
to the folder /Parameters of the model version, with a user defined identifier specified by the **model.ident** argument
in the **e2e_read()** function. The histories of proposed and accepted parameter combinations
are saved as csv files in the results folder and returned as a list object. The two csv files generated by the procedure have names:  
annealing_HRmult_proposalhistory-*, annealing_HRmult_acceptedhistory-*, where * denotes the value
of model.ident defined in the preceding e2e_read() function call.

\

To preserve the new harvest ratio multipliers and incorporate them into the fishing fleet model parameterisation
the multiplier values need to be applied to the scaling parameters which link the integrated effort by each gear to the harvest ratio
value which gets piped into the ecology model. Manually update the values in rows 12-21 (excluding the header row) of the
file /Parameters/fishing_fleet_parameters_*.csv, by multiplying the existing values by the new multipliers emerging from the annealing process. 
Then make sure to use a version of the multiplier file with all values set to 1.0 for all future runs.


\

If the edited file 'fishing_fleet_parameters_*.csv' is saved with a new identifier (*) then in order to use it in a subsequent run of the
StrathE2E model (using the e2e_run() function) it will be necessary to edit the MODEL_SETUP.csv file in the relevant /Models/variant folder to point to the new file.

\

The returned list object contains three dataframes: parameter_proposal_history, parameter_accepted_history, new_parameter_data (a list of three).
The proposal and accepted histories can be further analysed with the function e2e_plot_opt_diagnostics() to assess the performance of the optimization process.

\

```{r annealingcode2, eval = FALSE}
# Copy the 2003-2013 version of the North Sea model supplied with the package into a user workspace (Windows OS):
	e2e_copy("North_Sea", "2003-2013",
	dest.path="C:/Users/username/Documents/Models")
#
# Load the 2003-2013 version of the North Sea model from the user workspace:
	model<-e2e_read(model.name="North_Sea",
                          model.variant="2003-2013",
                          model.ident="TEST",
                          user.path="C:/Users/username/Documents/Models")
#
# Quick demonstration of the annealing function in operation:
e2e_optimize_hr(model, nyears=5, n_iter=10, start_temperature=0.5)
#
# More realistic configuration would be as follows :
e2e_optimize_hr(model, nyears=50, n_iter=1000, start_temperature=1)
# (WARNING - this will take about 26 hours to run)
#
```

\

\

## 5.6 Estimating of fishing activity rates **e2e_optimize_act()**

Optimize Strathe2E fishing gerar activity multipliers to maximize the likelihood of either observed ecosystem target data or known harvest ratios.

**Table 20**. Arguments of the function **e2e_optimize_act()**.

| Argument          | Description                                                                                                                                                                                                                               |
| ----------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| model             | R-list object generated by the e2e_read() function which defines the model configuration                                                                                                                                                  |
| selection         | Text string from a list to select the method for optimization. Select from: "HR", "ECO", Remember to include the phrase within "" quotes (default = "HR"). "HR" selects a scheme for optimizing the activity multipler to known harvest ratios and uses just the fishing fleet model. "ECO" selects a scheme for optimizing to the observational data in the state of the ecosystem and so uses both the fishing fleet and ecology models |
| n_iter            | Number of iterations of the model (default=500)                                                                                                                                                                                           |
| start_temperature | Initial value of the simulated annealing temperature parameter (default=0.5). Suggested values in the range 0.0005 - 5. Higher values increase the probability of rejecting parameter combinations producing an improvement in likelihood |
| cooling           | Rate at which the simulated annealing temperature declines with iterations (default=0.975). Suggested values in the range 0.9 - 0.985                                                                                                     |
| csv.output        | Logical. If FALSE then disable writing of csv output files - useful for testing (default=TRUE)                                                                                                                                            | 
| runtime.plot      | Logical. If FALSE then disable runtime plotting of the progress of the run - useful for testing (default=TRUE)                                                                                                                            |
| n_traj            | Specific to selection="HR" : Number of repeats (trajectories) of the simulated annealing process (default=100).                                                                                                                           |
| deltaHi           | Specific to selection="HR" : Initial coefficient of variation for jiggling the activity multiplier values (default=0.2). This is set in the parameter file "annealing_SD_fishing.csv" in the case where selection="ECO" and can be manually varied during a run |
| attenuationstep   | Specific to selection="HR" : Number of itereations between down-steps of the jiggling factor applied to the multiplier values. The jiggling rate is attenuated by a factor of 2 every xx iterations (default=500).                        |
| deltaGrat         | Specific to selection="HR" : Coefficient of variation for jiggling the gear linkage values (default=0.25).                                                                                                                                |
| nyears            | Specific to selection="ECO" : Number of years to run the ecology model in each iteration (default=40).                                                                                                                                    |
| quiet             | Logical. Specific to selection="ECO" : If TRUE then suppress informational messages at the start of each iteration (default=TRUE).                                                                                                        |

\

Usually, the fishing gear activity rates are a known input to the model. Sometimes, however, it may be necessary to use the model to hindcast the likely rates of 
activity given established values for the ecology model parameters, the scaling parameters linking effort to harvest ratios, and the relative spatial distribution of fioshing. 
This is a situation that might 
apply when the coupled fishing and ecology model system has been comprehensively fitted for a data-rich period of time, and is then being applied to a different 
data-poor period where the fishing gear activity rates are unknown or only approximately known.

\

Proposed values of the gear activity rates need to be contained in the standard **fishing_activity_parameters\*.csv** input file for the fishing fleet model. 
In addition, the model setup includes 
a separate csv file (**fishing_activity_multiplier.csv**) containing multipliers to be applied to each of these rates. These are useful when 
conducting scenario experiments, but would ordinarily be set to 1.0. In the same way as for the effort-harvest ratio scaling value fitting, the 
simulated annealing code for optimizing the gear activity rates iterates the multiplier values to seek the combination producing the best fit 
of the overall model to the set of target data. At the end of the process, the code saves the history of proposed and accepted multiplier values. 
If deemed reasonable, the final accepted version of the multipliers can then be manually applied to the corresponding rates in the initial input 
file to produce an updated version.

\

It is important to note that the **user.path** argument in the **e2e_read()** function call needs to point to a user workspace folder, not the default
North Sea model provided with the package. This is because the annealing function needs write-access to the model **/Parameters** folder,
but the **/extdata/Models** folder in the package installation is read-only. To use the annealing function on the North Sea model, use the **e2e_copy()** 
function to make a copy of the North Sea model in the user workspace.

\

There are two alternative ways to optimize fishing activity rates. The first is conceptually equivalent to the simulated annealing functions for
optimizing the ecology model parameters, in that the target for fitting is the same database of observational 
indices of ecosystem status. This uses the fully coupled fishing fleet and ecology model system and hence is costly to run.

\

The alternative function for estimating the activity rate of gears is much 
faster and relies only on the fishing fleet model, but here the target for optimization is a set of known harvest ratios. The circumstances under which 
the past harvest ratios for each guild are sufficiently well known to enable this fitting approach are likely to be rather restricted, but this method is 
probably the more robust. In both approaches the chief problem is that the selectivity patterns of the fishing gears are likely to be highly overlapping
so that it will be difficult to identify a realistic maximum likelihood set of activity rates. 

\

The target data for the second method for estimating activity rates are contained in the file **fishing_fleet_target_harvest_ratios\*.csv** (where * is a free text identifier) located in
the **/ModelName/VariantName/Target_data** folder of the model definition data.

\

The implementation of simulated annealing in this function is slightly different to the other optimization functions in the package. In this case,
because the problem of local maxima in the likelihood response surface is particularly acute (due to the potential overlap in the selectivities of different gears), so
the annealing process is replicated many times from different randomly selected initial conditions. We refer to these replicates as 'trajectories'. Each trajectory follows a 
different pathway through the parameter space. At the end of the process, the best-fit set of activity multipliers is selected from across all the trajectories. 
In addition, the coefficient for jiggling the parameters is systematically attenuated with increasing iterations rather than remaining constant or potentailly being manually
attenuated during the run my the user. Because the process uses only the fishing fleet model it is relatively fast so many trajectories can be run in a moderate time span.

\

To cope with the overlap in selectivity of the fishing gears, the function uses a set of linkage parameters specified in the the file **/Parameters/fishing_fleet_gear_linkages.csv**
These parameters force selected gears activities to vary in concert $\pm$ some variation, rather than independently. The parameters for the gear linkages are located in the file ../Modelname/Variantname/Parameters/fishing_fleet_gear_linkages.csv. 
The table of linkages specifies which gear activity rates are forced to 
vary in concert during the fitting process, as opposed to varying independently. The value of the linkage coefficient defines the 
scaling of changes in the activity rate of a dependent gear relative to its linked independent gear. For example, if gear 8 is permitted to vary 
independently (value in column "Gear to which linked" = NA and "Linkage coefficient" = NA). If gear 9 is dependent on gear 8 then the activity rate 
of gear 9 this would be specified by e.g. "Gear to which linked" = 8 and "Linkage coefficient" = 0.645. This would force the activity of gear 9 to be 
set at gear8_activity * (0.645 $\pm$ a random variation defined by the function argument **deltaGrat**).

\

The function produces a real-time graphical summary of the progress of the fitting procedure, displaying
the likelihoods of the proposed and accepted parameter sets at each iteration.
Y-axis (likelihood) range of the real time plot can be varied during the run by 
editing the setup file **annealing_SD_fishing.csv**

\

At the end of the procedure a new version of the gear activity multipliers file is exported
to the folder **/Parameters** of the model version, with a user defined identifier specified by the **model.ident** argument
in the **e2e_read()** function. The distributional statistics of the likelihood-weighted activity multipliers at the end of each trajectory, and the 
emergent harvest ratios relative to the target values are returned as a list object and optionally saved as csv files in the results folder (file names **fitting_to_knownHR_gearmult_stats-\*.csv** and 
**fitting_to_knownHR_HR_rel_target_stats\*.csv**, where * is the **model.ident** set by the **e2e_read()** function.

\

In order to use the new fitted parameter values in a subsequent run of the StrathE2E2 model (using the **e2e_run()** function)
it will be necessary to edit the **MODEL_SETUP.csv** file in the relevant **/Models/Variant** folder to point to the new harvest ratio multiplier file. Or, to 
preserve the new activity rate multipliers and incorporate them into the fishing fleet model parameterisation, manually update the values in the 
file **/Parameters/fishing_activity_parameters\*.csv**, by multiplying the existing values by the new multipliers emerging from 
the annealing process. If using Excel to do this task beware of the number formatting settings - specifically the number of decimal places displayed with scientific
notation. Increase the number of decimal places before saving as a csv file or precision of the activity rates will be lost.

\

If the edited file **fishing_activity_parameters\*.csv** is saved with a new identifier (*) then in order to use it in a subsequent run of the
StrathE2E2 model (using the **e2e_run()** function) it will be necessary to edit the **MODEL_SETUP.csv** file in the relevant **/Models/Variant** folder to 
point to the new file. Then make sure to use a version of the multiplier file with all values set to 1.0 for all future runs.

\

```{r annealingcode4, eval = FALSE}
# Copy the 1970-1999 version of the North Sea model supplied with the package into a user workspace (Windows OS):
e2e_copy("North_Sea", "1970-1999",
            dest.path="C:/Users/username/Documents/Models")
#
# Load the 1970-1999 version of the North Sea model from a user workspace:
model<-e2e_read(model.name="North_Sea",
                  model.variant="1970-1999",
                  model.ident="TEST",
                  user.path="C:/Users/username/Documents/Models")
# Quick Demonstration of optimizing to ecosystem data:
test_run <- e2e_optimize_act(model, selection="ECO", n_iter=10, start_temperature=0.4 , cooling=0.975, csv.output=FALSE, nyears=5 )
str(test_run,max.level=1)                         # View the structure of the returned list
#
# More realistic configuration would be (WARNING - this will take about 14 hours to run) :
find_gear_activity_mult_HR_target(model, n_traj=100, n_iter=3000, start_temperature=0.5)
test_run <- e2e_optimize_act(model, selection="ECO", n_iter=500, start_temperature=0.4 , cooling=0.975, csv.output=TRUE, nyears=40 )
# 
#
# Quick Demonstration of optimizing to harvest ratios:
test_run <- e2e_optimize_act(model, selection="HR", n_iter=100, start_temperature=1.0, cooling=0.985, csv.output=FALSE, n_traj=5 )
str(test_run,max.level=1)                         # View the structure of the returned list
str(test_run$new_parameter_data,max.level=1)      # View the structure of the returned list element containing parameter objects
test_run$new_parameter_data                       # View the new, final accepted parameter data#
#
# More realistic configuration would be (WARNING - this will take about 6 hours to run) :
model<-e2e_read(model.name="North_Sea",
                model.variant="1970-1999",
                model.ident="fittingrun",
                user.path="C:/Users/username/Documents/Models")
fitting_data<-e2e_optimize_act(model, selection="HR", n_iter=3000, start_temperature=1.0, cooling=0.985, csv.output=TRUE, n_traj=100 )
```

\

\

## 5.7 Diagnostics for optimization processes **e2e_plot_opt_diagnostics()**

Plot diagnostic data on the performance of parameter optimization functions.

**Table 21**. Arguments of the function **e2e_plot_opt_diagnostics()**.

| Argument          | Description                                                                                                                                                                                                                               |
| ----------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| model             | R-list object generated by the **e2e_read()** function which defines the model configuration                                                                                                                                                  |
| selection         | Text string from a list to select the method for optimization. Select from: "ECO", "HR", "ACT", corresponding to output generated by the e2e_optomize_eco(), e2e_optomize_hr(), or e2e_optomize_act() function, Remember to include the phrase within "" quotes (default = "HR"). "HR" selects a scheme for optimizing the activity multipler to known harvest ratios and uses just the fishing fleet model. "ECO" selects a scheme for optimizing to the observational data in the state of the ecosystem and so uses both the fishing fleet and ecology models |
| fitted.to         | Specific to the case where selection="ACT"; text string from a list identifying which version of activity optimization procedure generated the data to be plotted. Select from: "ECO", "HR", corresponding to maximising the likelihod of ecosystem state data, or zonal harvest ratios respectively. Default="". Remember to include the phrase within "" quotes |
| use.saved         | Logical. If TRUE use data from a prior user-defined run held as csv files data in the current results folder (default=FALSE).
| use.example       | Logical. If TRUE use pre-computed example data from the internal North Sea model rather than user-generated data (default=FALSE).
| results           | R-list object of output from an optimization process generated by the e2e_optomize_eco(), e2e_optomize_hr(), or e2e_optomize_act() function. Only needed if use.saved1=FALSE and use.example1=FALSE. (Default=NULL).

\
Creates diagnostic plots from outputs of the functions for optimizing ecology, harvest ratio, and fishing activity parameters: e2e_optimize_eco(), e2e_optimize_hr(), and e2e_optimize_act() respectively. 
 
\

The function takes the history of parameter values proposed during each simulated annealing process, expresses these relative to the initial value of each parameter, generates qunatiles (0.005, 0.25, 0.5, 0.75 and 0.995) of the distribution for each parameter, and plots the results as box-and-whisker diagrams.
The relativity to initial values is expressed as (proposed - initial)/initial so the initial is represented by a relative value of zero. For each parameter, the final accepted value is over-plotted onto the box-and-whisker as a red bar. The resulting diagram, with a separate box-and-whisker for each parameter
shows the extent to which each parameter has migrated from its initial value en-route to the final accepted state.

\

The format of the parameter diagram is the same for each type of optimization scheme: ecology, harvest ratio, and fishing activity parameters. However, in the case where activity parameters are
optimized to zonal harvest ratios rather than ecosystem state data, the diagram has an additional panel showing the distribution of proposal harvest ratios relative to the target, and (in red) the harvest 
ratios corresponding to the final accepted activity parmeters (relative to the targets).

\

Arguments for this function permit the input data to be drawn from an existing data object generated by the various optimization functions, previously genearted csv files, or example data
provided with the package for versions of the internal North Sea models.

\

Documentation in a dataframe format on each of the classes of parameters in the model can be obtained with the function **e2e_get_parmdoc()**. This provides a key to the abbreviated parameter names especially in the diagnostic plots for **e2e_optimize_eco()**.

\

As well as drawing the plot the function returns a list object containing a) an array of the quantiles of the proposal distributikn for each parameter, and b) an array of one row, of the final accepted parameter set from the procedure.

\

```{r optdiagnostics1, eval = FALSE}
# Load the 1970-1999 version of the North Sea model supplied with the package and generate a quick test data object
# with only 10 itereations and running the model for 5 years. More realistic would be at least 500 iterations and 
# running for 50 years. Even so this example will take a few minutes to run:
model<-e2e_read(model.name="North_Sea",
                model.variant="1970-1999",
                model.ident="test")
# This model is already optimized to the observed ecosystem data supplied with the package
# Perturb the temperature driving to knock the model away from its maximum likelihood state relative to the target data:
model$data$physics.drivers$so_temp <- model$data$physics.drivers$so_temp+3   # add 3 degC to upper layer offshore temperatures
model$data$physics.drivers$si_temp <- model$data$physics.drivers$si_temp+3   # add 3 degC to inshore temperatures
model$data$physics.drivers$d_temp  <- model$data$physics.drivers$d_temp+3    # add 3 degC to lower layer offshore temperatures
test_run  <- e2e_optimize_eco(model, nyears=5, n_iter=10, start_temperature=0.4,csv.output=FALSE,runtime.plot=FALSE)
plot_data <- e2e_plot_opt_diagnostics(model,selection="ECO",results=test_run)
str(plot_data,max.level=1)  # show the structure of the list object plot_data
```

```{r optdiagnostics2, fig.width=6, fig.height=8, fig.align='center'}
# Or... plot example date supplied with the package showing some data generated during the process of optimizing the North Sea model:
model     <- e2e_read(model.name="North_Sea", model.variant="1970-1999")
plot_data <- e2e_plot_opt_diagnostics(model,selection="ECO",use.example=TRUE)
```

**Figure 11**. Diagnostic plot of the procedure from optimizing ecology model parameters for the North Sea 1970-1999 model.

\

```{r optdiagnostics3, eval = FALSE}
# Same for harvest ratio optimization...
model<-e2e_read(model.name="North_Sea",
                model.variant="1970-1999",
                model.ident="test")
# This model is already optimized to the observed ecosystem data supplied with the package
# Perturb the temperature driving to knock the model away from its maximum likelihood state relative to the target data:
model$data$physics.drivers$so_temp <- model$data$physics.drivers$so_temp+3   # add 3 degC to upper layer offshore temperatures
model$data$physics.drivers$si_temp <- model$data$physics.drivers$si_temp+3   # add 3 degC to inshore temperatures
model$data$physics.drivers$d_temp  <- model$data$physics.drivers$d_temp+3    # add 3 degC to lower layer offshore temperatures
test_run  <- e2e_optimize_hr(model, nyears=5, n_iter=10, start_temperature=0.4,csv.output=FALSE,runtime.plot=FALSE)
plot_data <- e2e_plot_opt_diagnostics(model,selection="HR",results=test_run)
str(plot_data,max.level=1)  # show the structure of the list object plot_data
```

```{r optdiagnostics4, fig.width=6, fig.height=6, fig.align='center'}
# Or... plot example date supplied with the package showing some data generated during the process of optimizing the North Sea model:
model     <- e2e_read(model.name="North_Sea", model.variant="2003-2013")
plot_data <- e2e_plot_opt_diagnostics(model,selection="HR",use.example=TRUE)
```

**Figure 12**. Diagnostic plot of the procedure for optimizing harvest ratio multiplers parameters for the North Sea 2003-2013 model.

\

```{r optdiagnostics5, fig.width=6, fig.height=6, fig.align='center'}
# For activity rate optimization relative to ecosystem data:
model<-e2e_read(model.name="North_Sea",
                model.variant="1970-1999",
                model.ident="test")
# This model is already optimized to the observed ecosystem data supplied with the package, but not by optimizing gear activity rates
# The e2e_optimize_eco() function was used in this case.
# Perturb the temperature driving to knock the model away from its maximum likelihood state relative to the target data:
model$data$physics.drivers$so_temp <- model$data$physics.drivers$so_temp+3   # add 3 degC to upper layer offshore temperatures
model$data$physics.drivers$si_temp <- model$data$physics.drivers$si_temp+3   # add 3 degC to inshore temperatures
model$data$physics.drivers$d_temp  <- model$data$physics.drivers$d_temp+3    # add 3 degC to lower layer offshore temperatures
test_run  <- e2e_optimize_act(model, selection="ECO", n_iter=10, start_temperature=0.4, , cooling=0.975, csv.output=FALSE, nyears=5 )
plot_data <- e2e_plot_opt_diagnostics(model,selection="ACT",fitted.to="ECO",results=test_run)
str(plot_data,max.level=1)  # show the structure of the list object plot_data

# There are no example data available in the package for this function
```

\

```{r optdiagnostics6, fig.width=6, fig.height=6, fig.align='center'}
# For activity rate optimization relative to zonal harvest ratios:
model<-e2e_read(model.name="North_Sea",
                model.variant="1970-1999",
                model.ident="test")
# Activity rates in this model are already optimized to the target harvest ratios supplied with the package but we woud not expect 
# to recover these values in this short demonstration run
test_run  <- e2e_optimize_act(model, selection="HR", n_iter=100, start_temperature=1.0, , cooling=0.985, csv.output=FALSE, n_traj=5 )
plot_data <- e2e_plot_opt_diagnostics(model,selection="ACT",fitted.to="HR",results=test_run)
str(plot_data,max.level=1)  # show the structure of the list object plot_data
```

```{r optdiagnostics7, fig.width=6, fig.height=6, fig.align='center'}
# Or... plot example date supplied with the package showing some data generated during the process of optimizing the North Sea model:
model     <- e2e_read(model.name="North_Sea", model.variant="1970-1999")
plot_data <- e2e_plot_opt_diagnostics(model,selection="ACT",fitted.to="HR",use.example=TRUE)
```

**Figure 13**. Diagnostic plot of the procedure or optimizing gear activity multiplers parameters for the North Sea 1970-1999 model.

\

\

## 5.8 Calculating initial values for harvest ratio scaling parameters **e2e_calculate_hrscale()**

Calculate initial values of the scaling parameters for the fishing fleet model which link effort to harvest ratios.

**Table 22**. Arguments of the function **calculate_hr_scale_values()**.

| Argument          | Description                                                                                 |
| ----------------- | ------------------------------------------------------------------------------------------- |
| model             | R-list object generated by the e2e_read() function which defines the model configuration  |

\

The scaling parameters convert the effort applied by each gear to each living resource guild in the model, into a harvest ratio (proportion of mass captured per day).
Effort is the product of activity and the relative power of each gear with respect to each guild.

\

In order to estimate these scaling parameters, data are needed for a 'calibration' time period when activity, catch and harvest ratio are all known.

\

The function assumes that the relevant model configuration has already been loaded using the **e2e_read()** function. It is expected that
this has loaded the array of power parameters (catch per unit activity by gear and guild during the calibration period) given in the file /Parameters/fishing_power_parameters_*.csv,
the corresponding activity in the file /Parameters/fishing_activity_parameters_*.csv, and the known harvest ratios of each guild in the file /Target_data/known_regional_harvest_ratios_*.csv

\

The function returns rough estimates of the scaling parameters. These are only rough because they assume that the effort and resource biomass are distributed uniformly
over the model domain. Hence they should be used as initial estimates to be refined by optimization.

\

The scaling parameters are dispayed on the screen and saved in a dataframe. The values then need to be manually edited
into rows 12-21 (excluding the header row) of the file /Parameters/fishing_fleet_parameters_*.csv

\

```{r calculatecode1, eval = TRUE}
# Load the 2003-2013 version of the North Sea model supplied with the package and calculate scaling parameter values:
model <- e2e_read("North_Sea", "2003-2013")
#
scale_values <- e2e_calculate_hrscale(model)
scale_values
#
```

\

\

# 6 Sensitivity and Monte Carlo analysis

## 6.1 Background

### 6.1.1 Sensitivity analysis

Sensitivity analysis helps to highlight the parameters which have the most influence in the simulated annealing optimization process. 
In addition the analysis can help to identify parameters for which constraint by external information would be most beneficial.

\

Simple one-at-a-time sensitivity analysis of parameters in a model is a straightforward process when there are few parameters and 
relatively little interaction between parameters. However, StrathE2E2 include many parameters and there is a strong possibility of 
interactions, so a more sophisticated scheme is required.

\

The R package includes a function to perform global analysis of parameter sensitivity using the Morris method (Morris, 1991). This 
involves a factorial sampling scheme to generate parameter values for replicate model runs and generate values for the 'Elementary Effect' 
(EE) of each parameter. The mean and standard deviation of EE values over many runs allows ranking of the parameters in terms of the 
sensitivity of their effects on the model, and the strength of interactions with other parameters.

\

The parent parameter set for the sensitivity analysis ($\theta$) incorporates all of the parameters and inputs to the model, ideally. the 
maximum-likelihood set of fitted ecology model parameters produced by the simulated annealing scheme, plus the fixed ecological 
parameters, the physical configuration parameters of the model (layer thicknesses, inshore/offshore areas), the parameters passed 
to the ecology model from the fishing fleet model, and the fishing activity rates and environmental driving data. The sensitivity 
to the environmental driving data is investigated by applying a scaling factor uniformly to the annual cycle of environmental data values.

\

From this parent set, a series of child-sets ($\theta_{k}^{*}$ where $1 \ge k \ge r$) is generated by applying a separate random increment to each 
parameter; $\theta_{k}^{*} = \theta + \delta[k]$ where $\delta[k]$ is a vector of random values from a gaussian distribution of mean 0 and 
standard deviation 
given by a fixed coefficient of variation applied to the parent-set value of each parameter. This process is equivalent to that 
for generating parameter proposals in the simulated annealing scheme, except that all of the parameters and model inputs are varied 
in the sensitivity analysis, rather than just the subset of parameters to be fitted.

\

For each of the child-sets of parameters, the likelihood of the suite of observed data indices $P(observations|\theta_{k,0}^{*})$ is calculated 
following runs of the model to stationary state, exactly as in simulated annealing scheme. We refer to these as trajectories.

\

Then, for each trajectory, the parameters ($1 \ge i \ge n$) are varied in turn, one at a time, by adding a fixed proportionality increment, 
the model re-run, and the likelihood computed $(P(observations|\theta_{k,i}^{*}))$. The proportionality increment for a given trajectory is drawn 
at random from a set of four fixed levels in the range $\pm 10\%$ of the child parameter set (-10\%, -5\%, +5\%, +10\%, i.e. $\Delta$ = 0.9, 0.95, 1.05, 1.10), 
so that:

$$\theta[i]_{k,i}^{*} = \theta[i]_{k,0}^{*} . \Delta[z]$$

where $\theta[i]^{*}$ is the ith element of $\theta^{*}$, $\Delta[z]$ is the zth element of $\Delta$, and z is drawn at random from the series (1, 2, 3, 4) for each value 
of $k$. Hence, for each trajectory the model runs are repeated $n+1$ times, where i=0 corresponds to the baseline run for each trajectory. 
The total number of nested runs to support the sensitivity analysis is thus $r*(n+1)$.

\

From each level-run of the model (in which a single parameter is perturbed), and its corresponding likelihood, the Elementary Effect 
of the parameter is calculated as

$$EE_{k,i} = \frac{(P(observations|\theta_{k,i}^{*}) - P(observations|\theta_{k,0}^{*}))}{\sqrt{(\Delta[z]-1)^{2}}}$$


\

On completion of the runs for all trajectories, the mean ($\mu_{i}$) and standard deviation ($S_i$) of the r Elementary Effects for each parameter 
(i) are calculated. For the mean,

$$\mu_{i} = \frac{1}{r} \sum_{k=1}^{r} EE_{k,i}$$

\

The magnitude of the mean represents the sensitivity of each parameter, and the corresponding standard deviation indicates the degree of 
non-linearity in the response or interaction with other parameters. The standard error of the mean for each parameter ($SEM_i$) is given by:

$$SEM_i = \frac{S_i}{\sqrt{r}}$$

\

If $S_i > |\mu_i | . \frac{\sqrt{r}}{2}$ then we can approximately conclude that $\mu_i$ is significantly greater than zero.

\

\

### 6.1.2 Monte Carlo analysis

The Monte Carlo scheme provided with the StrathE2E2 package has two modes of operation referred to as
"baseline mode", and "scenario mode".

\

Baseline mode involves generating a list of ecology model parameter sets and, for each set, determining the
likelihood of observational data consistent with the environmental and fishery drivers. The model outputs
generated with each set are then weighted by the likelihood before deriving quantiles of their distribution.
The quantile ranges then represent credible intervals of the model outputs.

\

The ensemble of baseline mode parameter sets is generated by sampling the ecology model parameters from uniform distibutions
centred on an initial (baseline) parameter set. Ideally, this should be the maximum likelihood parameter set arrived at by prior application the various 
optimization functions in the package to 'fit' the model to a suite of observational data on the state of ecosystem
(target data for fitting are in /Modelname/Variantname/Target_data/annual_target_data_*.csv).
Each iteration in the Monte Carlo scheme then generates a unique time series of model outputs at daily intervals,
together with the overall likelihood of the obervational target data. From these data, the function then derives credible intervals of each output by
weighting the values from each parameter set by the associated likelihood.

\

For each of the new parameter sets, the likelihood of the suite of observed data indices $P(observation|\theta_{k}^{*})$ 
is calculated following a run of the 
model to stationary state (as described in section 5 on parameter optiization), and all the output from the final, stationary year of each run is 
saved ($k$ = 0 to $r$, where $k$ = 0 corresponds to the maximum likelihood (baseline) parameter set).

\

To calculate the credible interval for any direct or derived model output variable (e.g. annual average mass density of a state variable, or 
the mass density on a given day in the final year), the values from the individual model runs ($V_k$) and the associated likelihoods ($P_k$) are 
assembled as a list of ($r + 1$; $k$ = 0 to $r$) pairs ($a_k$ = ($V_k$,$P_k$)). The list for each variable is then sorted by ascending values $V$ 
($a_j$ ; $j = 1$ to ($r + 1$) such that $V_{j} \ge V_{(j-1)}$). Then, the cumulative likelihood 
and the cumulative likelihood with increasing variable value is calculated as follows

$$C_j = (P_{j=1}, \sum_{j=1}^{j=2} P_{j}, \sum_{j=1}^{j=3} P_{j}, \sum_{j=1}^{j=4} P_{j} ... \sum_{j=1}^{j=(r+1)} P_{j} )$$

Proportions of the maximum cumulative likelihood are then caclulated as $Q_{j} = C_{j} / C_{(r+1)}$

\

Finally, values of $V$ corresponding to discrete values of $Q$ = (0.005, 0.25, 0.5, 0.75, 0.995) are extracted by interpolation. These values 
span the 0.5\% and 99.5\% credible intervals of the model output given the observed target data and uncertainty in the fitted ecology parameters. 
Note that uncertainties in the fixed ecology parameters, fishing fleet parameters, and environmental driving data are not reflected in these credible intervals.

\

In the "scenario mode" of operation, there are no target data against which the outpiuts from each member of the ensemble can be compared. In this case, 
the parameter sets and their
associated likelihoods are taken from a pre-existing **baseline** mode simulation. The parameter sets are used to generate model outputs for a scenarios of
environmental or fishing drivers - e.g. increased activity by selected gears, or warmer sea temperatures.
These scenario outputs are then weighted by the **baseline** mode likelihoods before derivation of quantiles
and credible intervals.


\

\

## 6.2 Performing a sensitivity analysis **e2e_run_sens()**

Perform a parameter sensitivity analysis on the StrathE2E model.

**Table 23**. Arguments of the function **e2e_run_sens()**.

| Argument       | Description                                                                                                                                                                                                                                                                                                                                    |
| -------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| model          | R-list object generated by the e2e_read() function which defines the model configuration                                                                                                                                                                                                                                                     |
| nyears         | Number of years to run the model in each iteration (default=50)                                                                                                                                                                                                                                                                                |
| n_traj         | Number of trajectories of parameter sets (default=16)                                                                                                                                                                                                                                                                                          |
| trajsd         | Coefficient of variation used to set the standard deviation for the gaussian distribution from which new paramater values are drawn to create each trajectoy baseline from the initial parent values (default=0.0075)                                                                                                                          |
| n_setoflevels  | Number of fixed levels of coefficient of variation used to generate the individual parameter values in each level-run. Must be an even number (default=4)                                                                                                                                                                                      |
| v_setoflevels  | Maximum coefficient of variation for the set of levels (default=0.1, i.e. -10 percent to +10 percent)                                                                                                                                                                                                                                          |
| coldstart      | Logical. If TRUE then the run is starting from cold - which means that the first trajectory baseline is the parent configuration as specified in the 'model' list object. If FALSE then signifies that this is a parallel run whuch will later be merged with the 'coldstart=TRUE' run. In this case the first trajectory baseline is a derivative of the parent. Default=TRUE  |
| quiet          | Logical. If TRUE then suppress informational messages at the start of each iteration (default=TRUE)                                                                                                                                                                                                                                                   |
| postprocess    | Logical. If TRUE then process the results through to a final sorted list of parameter sensitivies for plotting. If FALSE just produce the raw results. The reason for NOT processing would be if the job has been shared across multiple machines/processors and several raw result files need to be merged before processing. Default=TRUE          |
| csv.output     | Logical. If FALSE then disable writing of csv output files - useful for testing (default=TRUE)          |
| runtime.plot   | Logical. If FALSE then disable runtime plotting of the progress of the run - useful for testing (default=TRUE)          |

\

Performs a one-at-a-time parameter sensitivity analysis on the StrathE2E model using
the Morris Method for factorial sampling of the physical configuration parameters, ecology model
parameters, the fishing fleet model parameters, and the
environmental forcings (see section 6.1.1. above for details).

\

The basis for the method is a scheme for sampling the model parameters, one at a time, from distributions around baseline sets, and testing the effects on the performance
of the model against some criterion. In this case, the criterion is the likelihood of the observational data on the state of the ecosystem that
are used as the target for parameter optimization in the various simulated annealing functions supplied with the package.

\

The process requires an initial set of parameters for the model. We refer to this as the 'parent' parameter set. It is recommended that this should be the parameters producing
the maximum likelihood of the observational target data (as estimated by e.g. the e2e_optimize_eco() function). The MODEL_SETUP.csv file in the folder
/Models/Modelname/Modelvariant/ should be configured to point to the relevant files, and then loaded with the e2e_read() function. 

\

From this parent set, a series of 'child' parameter sets are generated by applying a separate random increment to each parameter drawn from a
gaussian distribution of mean 0 and standard deviation given by a fixed coefficient
of variation applied to the parent-set value of each parameter. 

\

For each of the child-sets of parameters, the likelihood of the observed target data is calculated following runs of StrathE2E to stationary state.
We refer to these as trajectory baselines.

\

Then, for each trajectory, the parameters are varied in turn, one at a time, by adding a fixed proportionality increment to the trajectory baseline values, the model re-run,
and the likelihood computed. We refer to these as 'level runs'. The proportionality increment is the same for all of the level runs within a given trajectory, and is drawn at random from a set of fixed
levels distributed symetrically around 0 (e.g. -10, -5, +5, +10 percent, i.e. proportions of the trajectory baseline values = 0.9, 0.95, 1.05, 1.10).

\

For each level run, the 'Elementary Effect (EE)' of the given parameter is calculated from the difference between the level run likelihood and the corresponding trajectory baseline likelihood.

\

On completion of all the trajectories, the raw results are (optionally) post-processed to generate the mean and standard deviations of all the EE values for each parameter. EE_mean is an index of the magnitude
of the sensitivity, and EE_sd is an index of the extent of interaction with other parameters.

\

During the run the function produces a real-time plot for each trajectory, in which the x-axis represents the sequence of parameters, and the y-axis is the likelihood
of the target data. A horizontal red line indicates the likelihood of the parent parameter set, horizontal grey line indicates the likelihood for each trajectory baseline
and each level-run likelihood is shown by a symbol. The y-axis range can be changed in real-time by editing the setup file "/Models/Modelname/Modelvariant/Parameters/Parameter_SD_Control/OATsensitivity_SD.csv"

\

The outputs from the function are directed to csv files in the current "results" folder. The outputs are:
a) Table of parameter values applied in each run of the model (OAT_parameter_values-*.csv, where * = model.ident as defined by e2e_read())
b) Table of the likelihood and EE value for each trajectory/level run (OAT_results-*.csv)
c) If post-processing is selected, then a table of Mean EE and standard deviation of EE for each parameter, sorted by the absolute value of EE_mean (sorted_parameter_elemental_effects-*.csv)

\

WARNING - This function will take several days to run to completion on a single processor with even a modest number of iteration. 
The total number of model runs required to support the analysis is r*(n+1) where r is the number of trajectories and n is the number of parameters. 
The function incorporates all of the physical configuration parameters, fixed and fitted ecology model
parameters, the fishing fleet model parameters, and the environmental forcings into the analysis, so r = 453. Each model run needs to be sufficiently long to achieve a stationary state 
and as a consequnce a typical runtime of around 10h per trajectory. The mininimum recommended number of trajectories
is 15,  to the function can take several days to complete.

\

However, it is possible to spread the load over multiple processor/machines with 
arguments in the function allowing for management of this parallelisation. Afterwards, the raw results files are combined into a single data set
using the e2e_merge_sens_mc() function, and then processed using the function e2e_process_sens_mc(). 

\

A separate function e2e_plot_sens_mc() produces a graphical representation of the EE_mean and EE_sd results.

\

Example results data for both variants of the North Sea model generated with the default function settings are included in the package and can be plotted by
setting an argument in the **e2e_plot_sens()** function.

\

```{r sensitivityanalysis1, eval = FALSE}
# Load the 2003-2013 version of the North Sea model supplied with the package:
model <- e2e_read("North_Sea", "1970-1999")
#
# Run the sensitivity analysis process:
# WARNING - Running a full sensitivity analysis takes days of computer time on a single machine/processor because it involves a huge number of model runs.
# The example below is just a (relatively) quick minimalist demonstration and should NOT be taken as the basis for any analysis or conclusions.
# Even so, this minimalist demo could take 45 min to run to completion because it involves 1359 model runs.
sens_results <- e2e_run_sens(model, nyears=1, n_traj=3, csv.output=FALSE)
head(sens_results)
#
# A more realistic sensitivity analysis would be something like:
sens_results <- e2e_run_sens(model, nyears=50, n_traj=16, postprocess=TRUE)
# DO NOT launch this configuration unless you are prepared to wait many days for the results
#
# Example of parallelising the process:
# Launch two (or more) runs separately on different processors...
# Launch batch 1:
model1 <- e2e_read("North_Sea", "1970-1999", model.ident="BATCH1")
sens_results1 <- e2e_run_sens(model1, nyears=50, n_traj=10, coldstart=TRUE, postprocess=FALSE)   # On machine 1
#
model2 <- e2e_read("North_Sea", "1970-1999", model.ident="BATCH2")
sens_results2 <- e2e_run_sens(model2, nyears=50, n_traj=10, coldstart=FALSE, postprocess=FALSE)  # On machine 2
# Note that these two runs return only raw data since postprocess=FALSE
#
# Then, afterwards, merge the two raw results files with text-tags BATCH1 and BATCH2, and post process the combined file:
model3 <- e2e_read("North_Sea", "1970-1999", model.ident="COMBINED")
processed_data <- e2e_merge_sens_mc(model3, selection="SENS", ident.list<-c("BATCH1","BATCH2"), postprocess=TRUE, csv.output=TRUE)
#
# Plot a diagram of parameter sensitivities from the combined data
e2e_plot_sens_mc(model3, selection="SENS", use.example=FALSE)
```

\

\

## 6.3 Performing a Monte Carlo analysis **e2e_run_mc()**

Run a Monte Carlo simulation with StrathE2E and derive centiles of credible values of model outputs.

**Table 24**. Arguments of the function **e2e_run_mc()**.

| Argument              | Description                                                                                                                    |
| --------------------- | ------------------------------------------------------------------------------------------------------------------------------ |
| model                 | R-list object generated by the e2e_read() function which defines the model configuration                                       |
| nyears                | Number of years to run each instance of the model. Needs to be long enough to allow the model to attain a stationary state (default=50) |
| baseline.mode         | Logical. If TRUE then the simulation adopts a baseline mode. If FALSE then scenario mode (default=TRUE)                      |
| use.example.baseparms | Logical. Value required only if baseline.mode=FALSE. If TRUE then the baseline parameters set is drawn from example data provided with the package. If FALSE then a user generated baseline parameter set is expected |
| baseparms.ident       | Default = "". Value required if baseline.run=FALSE and use.example.baseparms=FALSE, in which case this is the model.ident string of the a user generated baseline parameter set. Remember to include the phrase within "" quotes |
| begin.sample          | Default = 1. Value required if baseline.mode=FALSE. Value sets the first row in the file of baseline parameter data from which the sequence of parameters sets to be used in this run will be taken. Set a value > 1 if this run is part of a parallel set of runs |
| n_iter                | Number of iterations of the model (default=1000). If baseline.mode=FALSE and the number of sets in the supplied baseline parameter file beyond the value of begin.sample is less than n_iter, then n_iter is reset to the remaining number available |
| csv.output            | Logical. If FALSE then disable writing of csv output files - useful for testing (default=TRUE)                                 |
| runtime.plot          | Logical. If FALSE then disable runtime plotting of the progress of the run - useful for testing (default=TRUE)                 |
| postprocess           | Logical. If FALSE then disable post-processing of the cumulative outputs, for example in parallel mode where multiple runs are going to be combined later (default=TRUE) |


\

The Monte Carlo scheme provided with the StrathE2E2 package has two modes of operation. The first is referred to as
"baseline mode". This involves generating a list of ecology model parameter sets and, for each set, determining the
likelihood of observational data consistent with the environmental and fishery drivers. The model outputs
generated with each set are then weighted by the likelihood before deriving quantiles of their distribution.
The quantile ranges then represent credible intervals of the model outputs.

\

The second mode of operation is referred to as the "scenario mode". In this case, the parameter sets and their
associated likelihoods from a baseline mode simulation, are used to generate model outputs for scenarios of
environmental or fishing drivers - e.g. increased activity by selected gears, or warmer sea temperatures.
These scenario outputs are then weighted by the baseline mode likelihoods before derivation of quantiles
and credible intervals.

\

In baseline mode, the Monte Carlo scheme generates an ensemble of model runs by sampling the ecology model parameters from uniform distibutions
centred on an initial (baseline) parameter set. Ideally, this should be the maximum likelihood parameter set arrived at by prior application the various 
optimization functions in the package to 'fit' the model to a suite of observational data on the state of ecosystem
(target data for fitting are in /Modelname/Variantname/Target_data/annual_target_data_*.csv).
Each iteration in the Monte Carlo scheme then generates a unique time series of model outputs at daily intervals,
together with the overall likelihood of the obervational target data. From these data, the function then derives credible intervals of each output by
weighting the values from each parameter set by the associated likelihood.

\

Running a scenario mode simulation requires a baseline model to have been previously completed, in order to generate a list of parameter sets and associated
likelihoods. In scenario mode, model driving data (e.g. temperatures or fishing activities) are varied from the baseline configuration by editing
the model object generated by a e2e_read() function call, and the model is then run with the existing baseline model parameter sets, and the results
weighted by the baseline mode likelihoods in order to derive the credible intervals.

\

The choice of baseline vs scenario mode of simulations is determined in the function by a logical argument setting 'baseline.mode'.
 
\

In baseline mode, the coefficients of variation for jiggling the ecology parameter can be varied in real-time
during the run by editing the file "CredIntSim_SD.csv" in the folder
/Parameters/Parameter_SD_control/ of the model version. However, it is recommended to leave the setting constant for the duration of a run.
A CV of 0.10 to 0.15 is recommended. If comparing the credible intervals for two different baseline models then it is important to use the same CV in both cases.

\

In scenario model, the parameter sets are pre-ordained by a prior baseline simulation, so the CV setting is ineffective.

\

On completion of all the iterations of the model, whether in baseline or scenarion mode, for each model output variable in turn, values from the individual runs and
their associated model likelihoods are assembled as a list of paired values. The list is then sorted by ascending values of the model variable,
and the cumulative likelihoods with increasing value of model variable is calculated. Values for the model variable at standard proportions
(0.005, 0.25, 0.5, 0.75, 0.995) of the maximum cumulative likelhood are then extracted by interpolation.
These represent the credible intervals of model output given uncertainty in the ecology model parameters.

\

Outputs from the e2e_run_mc() function depend on whether post-processing is selected. If not, then raw data are saved as csv files
(by default) and the function returns a dataframe of the parameter sets used in the run and their corresponding likelihoods. In scenario mode this is just a replica of the 
baseline mode data that have been used to run teh simulation. If post-processing is selected, then 
both the raw and processed credible interval data are saved to csv files (by dafault), and the function returns a list opject which includes the parameter
dataframe plus all the processed credible interval data structures. csv files are saved in the folder /results/Modelname/Variantname/CredInt/ with an
identifier for the simulation (model.ident) created by the e2e_read() function. There are two types of output. First is simply an
accumulation of all the standard outputs from StrathE2E on a run-by-run basis. The second is a set of files containing the derived credible intervals.

\

The function displays a real-time graphic to show the progress of the simulation. During the StrathE2E-running phase of the
process an x-y graph is updated after each iteration (x-axis=iterations, y-axis=Likelihood of the target data), with a horizontal grey line showing the
baseline (maximum liklihood) result, and black symbols showing the likelihood for each iteration based on the parameter values sampled from the baseline.
The y-axis limits for the graph can be varied in real-time
during the run by editing the file "CredIntSim_SD.csv" in the folder /Parameters/Parameter_SD_control/ of the model version.

\

During the post-processing phase, a message is displayed as each output variable is processed.

\

WARNING - the scheme can take a long time to run (~2 days with the default settings), and generate large output files (11 files total 3.2 Mb per iteration). 
Check the memory capacity of you machine before starting a long run. It can make sense to 
parallelise the process by splitting across different processors and merging the results afterwards. The package contains a function for merging cumulative output files from 
multiple runs and post-processing the combined data.

\

In baseline mode, parallel instances of the function can be implemented on different processors, all starting from the same
initial parameter set. During the subsequent merging process, the outputs from the first parameter set, of all but the first instance, are stripped away and discarded.

\

Impementing parallel runs in scenario mode requires a different approach. Here, the parameter sets are pre-ordained and it is important to avoid using duplicate sets in the simulations.
Hence, the function includes an argument 'begin.sample' to select a pointer in the baseline parameter set to begin sampling in any scenario mode run. For example, in the first instance, begin.sample=1, and the number
of iterations (n_iter) might be set to e.g. 200. For the second instance the pointer (begin.sample) would then be set to 201, and so on. If begin.sample > 1, the function runs the first baseline parameter set (sample 1)
first before proceeding to the pointer begin.sample and completing the requested number of iterations, so there will be one extra iteration in the output files. If the requested number of iterations 
means that sampling would over the end of available number of parameter sets
then the number of iterations is reduced accordingly.

\

Although the e2e_run_mc() function generates large raw data files (11 files total 3.2 Mb per iteration), the processed data are much smaller (13 files total ~4 Mb regardless of
the nunber of iterations).



```{r montecarloanalysis1, eval = FALSE}
# Load the 1970-1999 version of the North Sea model supplied with the package:
   model <- e2e_read("North_Sea", "1970-1999")
# Run the Monte Carlo process and generate some data but disable csv output
# A quick demonstration of baseline mode with postprocessing and csv output disabled:
   demo <- e2e_run_mc(model,nyears=2,baseline.mode=TRUE,n_iter=10,csv.output=FALSE)
      str(demo,max.level=1)                     # View the structure of the returned list
      str(demo$CI_annual_avmass,max.level=1)    # View the structure of a returned list element containing sub-sets
      str(demo$CI_annual_fluxes,max.level=1)    # View the structure of a returned list element containing sub-sets
      head(demo$CI_annual_avmass$whole)         # View the top few rows on the whole domain data on annual averge mass
#
# --------------------------------------------------------------------------------------
#
# A more meaningful run would be:
# WARNING: This run will take about 48h to complete, much better to split up and spread across multiple processors !
   mc_results <- e2e_run_mc(model,baseline.mode=TRUE,nyears=50,n_iter=1000)
#
# --------------------------------------------------------------------------------------
#
# A quick demonstration run in baseline mode and save the data to csv files:
   basemodel <- e2e_read("North_Sea", "1970-1999",model.ident="mcbaseline")
   basedemo <- e2e_run_mc(basemodel,nyears=2,baseline.mode=TRUE,
                          n_iter=10,csv.output=TRUE)
#
# --------------------------------------------------------------------------------------
#
# Then a quick demonstration run in scenario mode using the saved baseline parameter data, and save to csv:
# First create an extreme fishing scenario - quadruple some gear activities, run for 10 year
   scenariomodel<-basemodel
   scenariomodel$setup$model.ident <- "mcscenario"
   scenariomodel$data$fleet.model$gear_mult[1] <- 4 # Gear 1 (Pelagic trawls) activity rate rescaled to 4*baseline
   scenariomodel$data$fleet.model$gear_mult[4] <- 4 # Gear 4 (Beam_Trawl_BT1+BT2) activity rate rescaled to 4*baseline
   scendemo <- e2e_run_mc(scenariomodel,nyears=10,baseline.mode=FALSE, 
                          use.example.baseparms=FALSE, baseparms.ident="mcbaseline",begin.sample=1,
                          n_iter=10,csv.output=TRUE)
#
# Compare the results of the baseline and scenario:
   basemodel <- e2e_read("North_Sea", "1970-1999", model.ident="mcbaseline")
   scenariomodel <- e2e_read("North_Sea", "1970-1999",model.ident="mcscenario")
   e2e_compare_runs_box(selection="ANNUAL", model1=basemodel, ci.data1=TRUE, use.saved1=TRUE,
#                        model2=scenariomodel, ci.data2=TRUE, use.saved2=TRUE )
   dev.new()
   e2e_compare_runs_box(selection="MONTHLY", model1=basemodel, ci.data1=TRUE, use.saved1=TRUE,
                        model2=scenariomodel, ci.data2=TRUE, use.saved2=TRUE )
#
# --------------------------------------------------------------------------------------
#
# Example of parallelising the process in baseline mode:
# Launch batch 1:
   model1 <- e2e_read("North_Sea", "1970-1999", model.ident="BATCH1")     ### On processor 1
   results1 <- e2e_run_mc(model1,nyears=2,baseline.mode=TRUE,
                          n_iter=10,csv.output=TRUE,postprocess=FALSE)                
# Launch batch 2 (on a different processor):
   model2 <- e2e_read("North_Sea", "1970-1999", model.ident="BATCH2")     ### On processor 2
   results2 <- e2e_run_mc(model2,nyears=2,baseline.mode=TRUE,
                          n_iter=10,csv.output=TRUE,postprocess=FALSE) 
# Note that these two runs return only raw data since postprocess=FALSE
# Note that Batch 2 requests 11 iterations, rather that 10 in Batch 1. The number of iterations do not
# have to be the same in each batch. However, the first in each batch has to use the initial parameter set 
# from the model setup, as this is the parent for all the subsequent samples generated during the run. When we
# come to merge the data from separate batches, data from the first iteration are stripped off and discarded for
# all but the first batch as we do not want to include duplicate data in the combined files. Hence we choose 11
# iterations here in Batch 2 to make the point, and we expect the combined data to include 20 iterations.

# Then, afterwards, merge the two raw results files with text-tags BATCH1 and BATCH2, and post process the combined file:
   model3 <- e2e_read("North_Sea", "1970-1999", model.ident="COMBINED")
   processed_data <- e2e_merge_sens_mc(model3, selection="MC", ident.list<-c("BATCH1","BATCH2"), postprocess=TRUE, csv.output=TRUE)
# or...
   combined_data <- e2e_merge_sens_mc(model3, selection="MC", ident.list<-c("BATCH1","BATCH2"), postprocess=FALSE, csv.output=TRUE)
   processed_data<-e2e_process_sens_mc(model3, selection="MC", use.example=FALSE,csv.output=TRUE)
#
# Plot the parameter likelihood distributions from the combined data
   e2e_plot_sens_mc(model3, selection="MC")
#
# --------------------------------------------------------------------------------------
#
# Example of parallelising the process in scenario mode, using the baseline parameter sets 'COMBINED' from above:
# The activity of all fishing gears is reduced to zero to create a no-fishing scenario.
# Run each batch for 10 years as a relatively quick demo - a real run would need to run for at least 40 year
# Launch batch 1:
   model1s <- e2e_read("North_Sea", "1970-1999", model.ident="BATCH1_S")     ### On processor 1
   model1s$data$fleet.model$gear_mult[1:12] <- 0 # Activity rates of all 12 gears rescaled to 0*baseline
   results1s <- e2e_run_mc(model1s,nyears=10,baseline.mode=FALSE, baseparms.ident="COMBINED", begin.sample=1,
                          n_iter=10,csv.output=TRUE,postprocess=FALSE)                
# Launch batch 2 (on a different processor):
   model2s <- e2e_read("North_Sea", "1970-1999", model.ident="BATCH2_S")     ### On processor 2
   model2s$data$fleet.model$gear_mult[1:12] <- 0 # Activity rates of all 12 gears rescaled to 0*baseline
   results2s <- e2e_run_mc(model2s,nyears=10,baseline.mode=FALSE,, baseparms.ident="COMBINED", begin.sample=11,
                          n_iter=10,csv.output=TRUE,postprocess=FALSE)
# Note that Batch 1 samples rows 1:10 of the baseline mode parameter set archive "COMBINED" (begin.sample=1, n_iter=10),
# so Batch 2 needs to start sampling at row 11 (begin.sample=11). The baseline archive contains 20 rows, so Batch 2 
# has the capacity to undertake up to 10 sample iterations (rows 11:20). If we select more than 10 iterations 
# (e.g. n_iter=15) then the code will automatically restrict to 10. Note that in fact, to be consistent with the 
# format of output files from the baseline mode, each scenario mode run where 'begin.sample' > 1 will complete 
# n_iter+1 iterations, by adding an initial run using the parameter values from row 1 of the baseline parameter
# set - which is then stripped off during merging. 
#
# Then, merge the two raw results files with text-tags BATCH1_S and BATCH2_S, and post process the combined file:
   model3s <- e2e_read("North_Sea", "1970-1999", model.ident="COMBINED_S")
   processed_datas <- e2e_merge_sens_mc(model3s, selection="MC", ident.list<-c("BATCH1_S","BATCH2_S"), postprocess=TRUE, csv.output=TRUE)
#
# Finally plot comparisons of the baseline and scenario model runs:
   e2e_compare_runs_box(selection="ANNUAL", model1=model3, ci.data1=TRUE, use.saved1=TRUE,
                        model2=model3s, ci.data2=TRUE, use.saved2=TRUE )
   dev.new()
   e2e_compare_runs_box(selection="MONTHLY", model1=model3, ci.data1=TRUE, use.saved1=TRUE,
                        model2=model3s, ci.data2=TRUE, use.saved2=TRUE )
#
```

\

\

## 6.4 Parallelisation of sensitivity and Monte Carlo analyses

Both the sensitivity and Monte Carlo analyses are extremely time consuming, especailly the sensitivity analysis which can take several day to 
complete a meaningful number of trajectories. This is beacause of the 
large number of replicate simulations that are required. However, the process is straightforward to parallelise on on multiple processors 
with the results being merged afterwards for processing.

\

### 6.4.1 Merging data generated on parallel processors **e2e_merge_sens_mc()**

Combine two or more sets of raw output data from parallel sensitivity or Monte Carlo analysis runs performed on separate machines/processors.

**Table 25**. Arguments of the function **e2e_merges_sens_mc()**.

| Argument        | Description                                                                                                                  |
| --------------- | ---------------------------------------------------------------------------------------------------------------------------- |
| model           | R-list object generated by the **e2e_read()** function which defines the model configuration                                 |
| selection       | Text string from a list to select the method for optimization. Select from: "SENS" or "MC", corresponding to output generated by the e2e_optomize_eco(), e2e_optomize_hr(), or e2e_optomize_act() function, Remember to include the phrase within "" quotes (default = "HR"). "HR" selects a scheme for optimizing the activity multipler to known harvest ratios and uses just the fishing fleet model. "ECO" selects a scheme for optimizing to the observational data in the state of the ecosystem and so uses both the fishing fleet and ecology models |
| ident.list      | A vector of text variables corresponding to the "model.ident" identifiers for each of the files to merged (list must be length 2 or greater) |
| postprocess     | Logical. if TRUE then process the results through to final data; if FALSE just produce the combined raw results (default=TRUE). The reason for NOT processing would be if there are further run results stil to be combined with the set produced by this function |
| csv.output      | Logical. If FALSE then disable writing of CSV output files - useful for testing (default=TRUE) |

\

The process of merging results files generated oon parallel processors is not as simple as merely concatenating the files as it is necessary to keep track of 
the unique identities of the sets of trajectories.

\

The files to be combined must be transferred into the same folder, and this is where the new combined files will be placed.
The path to locate the files is set in a e2e_read() function call. If not specified it is assumed that the files are located in the
default /results/Modelname/Variantname folder in the current user workspace.

\

An identifying text string for the new combined files is set by the 'model.ident' argument in a e2e_read() function call.

\

The list of files to be combined (any number > 1) is defined by a vector of their individual "model.ident" identifiers ("ident.list" argument). 

\

When combining the files, the function creates a seamless sequence of trajectory identifiers through the combined data, beginning from 1 for the first baseline (maximum likelihood) trajectory
of the first set.

\

If for any reason there is a need to combine separate batches of multiple run results, then 
post-processing can be delayed with the 'postprocess' argument until the last merge when all the data have been gathered
together. Stand-alone postprocess can be performed using the function e2e_process_sens_mc().

\

The function returns csv files of the merged data in with thw identifier model.ident as specified in a prior call of the e2e_real() function. 
If the argument postprocess=TRUE then a dataframe of processed output (for sensitivity analysis) or a list object of dataframes (for Monte Carlo analysis) 
is also returned.

\

**Details relating to merging sensitivity analysis files:**

e2e_run_sens() generates two output files per run - OAT_results-*.csv, and OAT_parameter_values-*.csv, where * is
a model.ident text string set as an argument of the e2e_read() function. Thus function merges both of these types of files.

\

When combining sensitivity analysis data the first-named model.ident in the ident.list vector MUST correspond to a run of the e2e_run_sens() function with the argument coldstart=TRUE, and all others with coldstart=FALSE. This forces the
first trajectory of sensitivity test to be performed on the baseline (e.g. maximum-likelihood) parameter set loaded with the initial e2e_read() function call. Thus is important
for the post-processing stage of the analysis which needs to be performed on the combined results.

\

**Details relating to merging Monte Carlo analysis files:**

e2e_run_mc() generates 11 different output files of accumulated outputs from the iterations (total 3.2 Mb per iteration) covering the range of model outputs. This merging function combines batches of each of these types of files. 
Check the memory capacity of you machine before starting a long run or before merging runs. The processed output consists of 13 files total ~4 Mb regardless of the number of iterations.

\

```{r sensmcmerge1, eval = FALSE}
# Example of parallelising the sensitivity analysis process:
# Launch two (or more) runs separately on different processors...
# Launch batch 1:
   model1 <- e2e_read("North_Sea", "1970-1999", model.ident="BATCH1")
   sens_results <- e2e_run_sens(model1, nyears=50, n_traj=10, coldstart=TRUE, postprocess=FALSE)   # On processor 1
# Note that coldstart=TRUE for the first batch only.
# Launch batch 2 (on a different processor):
   model2 <- e2e_read("North_Sea", "1970-1999", model.ident="BATCH2")
   sens_results <- e2e_run_sens(model2, nyears=50, n_traj=10, coldstart=FALSE, postprocess=FALSE)  # On processor 2
# Note that these two runs return only raw data since postprocess=FALSE

# Then, afterwards, merge the two raw results files with text-tags BATCH1 and BATCH2, and post process the combined file:
   model3 <- e2e_read("North_Sea", "1970-1999", model.ident="COMBINED")
   processed_data <- e2e_merge_sens_mc(model3, selection="SENS", ident.list<-c("BATCH1","BATCH2"), postprocess=TRUE, csv.output=TRUE)
# or...
   combined_data <- e2e_merge_sens_mc(model3, selection="SENS", ident.list<-c("BATCH1","BATCH2"), postprocess=FALSE, csv.output=TRUE)
   processed_data <- e2e_process_sens_mc(model3, selection="SENS", use.example=FALSE,csv.output=TRUE)
#
#
# Example of parallelising the Monte Carlo process:
# Launch batch 1:
   model1 <- e2e_read("North_Sea", "1970-1999", model.ident="BATCH1")     ### On processor 1
   results1 <- e2e_run_mc(model1,nyears=2,baseline.mode=TRUE,
                          n_iter=10,csv.output=TRUE,postprocess=FALSE)                
# Launch batch 2 (on a different processor):
   model2 <- e2e_read("North_Sea", "1970-1999", model.ident="BATCH2")     ### On processor 2
   results2 <- e2e_run_mc(model2,nyears=2,baseline.mode=TRUE,
                          n_iter=11,csv.output=TRUE,postprocess=FALSE) 
# Note that these two runs return only raw data since postprocess=FALSE
# Note 11 iterations in batch 2 - the first iteration wil be stripped off at merging so the
# combined data should include only 20 iterations.

# Then, afterwards, merge the two raw results files with text-tags BATCH1 and BATCH2, and post process the combined file:
   model3 <- e2e_read("North_Sea", "1970-1999", model.ident="COMBINED")
   processed_data <- e2e_merge_sens_mc(model3, selection="MC", ident.list<-c("BATCH1","BATCH2"), postprocess=TRUE, csv.output=TRUE)
# or...
   combined_data <- e2e_merge_sens_mc(model3, selection="MC", ident.list<-c("BATCH1","BATCH2"), postprocess=FALSE, csv.output=TRUE)
   processed_data<-e2e_process_sens_mc(model3, selection="MC", use.example=FALSE,csv.output=TRUE)

#
```

\

\

### 6.4.2 Processing data generated on parallel processors **e2e_process_sens_mc()**

Post-process already-created data from sensitivity or Monte Carlo analysis runs.

**Table 26**. Arguments of the function **e2e_merges_sens_mc()**.

| Argument        | Description                                                                                                                  |
| --------------- | ---------------------------------------------------------------------------------------------------------------------------- |
| model           | R-list object generated by the **e2e_read()** function which defines the model configuration                                 |
| selection       | Text string from a list to select the method for optimization. Select from: "SENS" or "MC", corresponding to output generated by the e2e_optomize_eco(), e2e_optomize_hr(), or e2e_optomize_act() function, Remember to include the phrase within "" quotes (default = "HR"). "HR" selects a scheme for optimizing the activity multipler to known harvest ratios and uses just the fishing fleet model. "ECO" selects a scheme for optimizing to the observational data in the state of the ecosystem and so uses both the fishing fleet and ecology models |
| ident.list      | A vector of text variables corresponding to the "model.ident" identifiers for each of the files to merged (list must be length 2 or greater) |
| use.example     | Logical. If TRUE use pre-computed example data from the internal North Sea model rather than user-generated data (default=FALSE) |
| csv.output      | Logical. If FALSE then disable writing of CSV output files - useful for testing (default=TRUE) |

\

The functions **e2e_run_sens()** and **e2e_run_mc()** are extremely time consuming so it makes sense to share the 
load across multiple processors in parallel and combine the results afterwards using the function **e2e_merge_sens_mc()**. Both the original
run functions and the merging function have post-processing options in their workflow. Nevertheless, there may be occasions where it is
necessary to simply post-process an already created raw data set. Hence the need for this function.

\

The model.ident identifier and folder path for input files must be set in a **e2e_read()** function call. If enabled, output files will be directed to the same folder.

\

**Details relating to sensitivity analysis processing:**

The function reads a file of raw sensitivity analysis data (OAT_results-*.csv, where * refers to the identifier model.ident set in the **e2e_read()** function) from a /results folder
and performs the post-processing that would ordinarily be done automatically within the e2e_run_sens() or e2e_merge_sens_mc() functions. 
The output of the processing is a dataframe (and optionaly a csv file) containing the mean Elementary Effect (EE_mean) and the standard deviation of EE (EE_sd) for each parameter, sorted by the absolute value of EE_mean.
EE_mean is an index of the magnitude of the sensitivity of a parameter, and EE_sd is an index of the extent of interaction with other parameters.

\

csv output from the function is a file named sorted_parameter_elementary_effects-*.csv in the current folder /results/Modelname/Variantname/

\

Optionally, the function can read example raw sensitivity data for one of the two North Sea model variants supplied with the package.

\

For details of how the Elementary Effect values are derived for each parameter see help(e2e_run_sens)

\

**Details relating to Monte Carlo analysis processing:**

The function **e2e_run_mc()** generates 11 different output files covering the range of data generated by StrathE2E. The post processing reads all of these and generates quantiles of the 
likelihood weighted distributions of all of the outputs and saves these as credible interval files.

\

By default, the function writes output to csv files, but this can be disabled by a function argument. The function always returns all of the processed data as
a list object. The function takes a few minutes to run as it has a lot of work to do, especially in processing all of the daily output variables.

\

Although the e2e_run_mc() function generates large raw data files (11 files total 3.2 Mb per iteration), the processed data are much smaller (13 files total ~4 Mb regardless of
the nunber of iterations.


```{r sensmcprocess1, eval = FALSE}
# Sensitivity analysis processing:
# Load details of the 1970-1999 version of the North Sea model supplied with the package:
model <- e2e_read("North_Sea", "1970-1999")
# Process the example data for this model variant provided with the package
   sens_results <- e2e_process_sens_mc(model, selection="SENS", use.example=TRUE)
   head(sens_results)
#'
# Monte Carlo analysis processing
# The example here assumes that raw results data have been previously
# generated by a of the e2e_run_mc() function, with the model.ident value
# assigned as "testdata", or that data from parallel runs have been gathered
# together in the same folder with model.ident="testdata".
# Load the model name and version for the data to be processed :
# e.g. ... model <- read_model("North_Sea", "2003-2013", model.ident="testdata")
   mc_results <- e2e_process_sens_mc(model, selection="MC", csv.output=TRUE)
   str(mc_results,max.level=1)
#
```

\

\

## 6.5 Diagnostics for sensitivity and Monte Carlo analyses **e2e_plot_sens_mc()**

Plot diagnostic data from sensitivity and Monte Carlo analyses.

**Table 27**. Arguments of the function **e2e_plot_sens_mc()**.

| Argument        | Description                                                                                                                  |
| --------------- | ---------------------------------------------------------------------------------------------------------------------------- |
| model           | R-list object generated by the **e2e_read()** function which defines the model configuration                                 |
| selection       | Text string from a list to select the method for optimization. Select from: "SENS" or "MC", corresponding to output generated by the e2e_optomize_eco(), e2e_optomize_hr(), or e2e_optomize_act() function, Remember to include the phrase within "" quotes (default = "HR"). "HR" selects a scheme for optimizing the activity multipler to known harvest ratios and uses just the fishing fleet model. "ECO" selects a scheme for optimizing to the observational data in the state of the ecosystem and so uses both the fishing fleet and ecology models |
| use.example     | Logical. If TRUE use pre-computed example data from the internal North Sea model rather than user-generated data (default=FALSE) |

\

Generates diagnostic plots from either a parameter sensitivity analysis or a Monte Carlo simulation of credible intervals of model outputs.

\

For both sensitivity and Monte Carlo analyses the inputs required by the function are a model list object created by the **e2e_read()** function defining 
the model configuration and relevant model.ident
argument to point to the required files, and a logical argument to determne the origin of data files (either saved as csv files in the user workspace, or
example data for the North Sea model provided with the package.

\

**Details relating to sensitivity analysis plot:**

The function reads processed results generated by the function **e2e_run_sens()**, or the function **e2e_process_sens_mc(...,selection="SENS")**, and plots the 
Elementary Effect mean (x-axis; magnitude of sensitivity) against Elementary Effect standard deviation (y-axis; strength of interactions between parameters).

\

Processed results from the function **e2e_run_sens()**, or the function **e2e_process_sens_mc(...,selection="SENS")**, are stored in the csv file
/results/Modelname/Variantname/sorted_parameter_elementary_effects-*.csv where * represents the model run identifier (model.ident) text embedded in the 
R-list object created by the e2e_read() function.

\

Each symbol in the plot represents a single parameter in the model. The parameters are colour-coded to indicate 6 different types - fitted and fixed
parameters of the ecology model, fishing fleet model parameters, fishery harvest ratios, environmental drivers, and physical configuration parameters.

\

The plot also shows a wedge formed by the two dashed lines. These correspond to $\pm$ 2 standard errors of the mean, so for points falling outside of the
wedge there is a significant expectation that the distribution of elementary effects is non-zero. For points falling within the wedge the distribution of
elementary effects is not significantly different from zero.

For details of how the Elementary Effect values are derived for each parameter see help(e2e_run_sens)

\

**Details relating to Monte Carlo analysis plot:**

The function creates a plot showing the credible distributions of ecology model parameters based on the results from the e2e_run_mc() function.
These distributions are formed from the input distributions to the Monte Carlo process, weighted by the likelihood of observed target data
on the state of the ecosystem given each combination of parameter values.

\

Post-processed data from the **e2e_run_mc()** function are stored in the file
/results/Modelname/Variantname/CredInt/CredInt_processed_PARAMETERresults_data-*.csv, where * represents the model run identifier
(model.ident) text embedded in the R-list object created by the e2e_read() function.

\

Each parameter in the plot is scaled to its baseline value as (value - baseline)/baseline. Ideally, the baseline is the maximum likelihood model
developed by application of the optimization functions available in the package (e.g. e2e_optimize_eco()).
Each parameter is then represented by a box and whisker plot which shows the distribution of credible parameter values around zero, i.e.
around the baseline. The median of the credible values distribution for each parameter
is shown my a black tick-mark. The box spans the 50% credible interval (quartiles of the cumulative likelihood
of simulated values). Whisker lines span the 99% credible interval.

\

The individual parameters are identified by numbers (rather than text names). These numbers correspond to the column numnbers in
the file /results/Modelname/Variantname/CredInt/CredInt_processed_parameters-*.csv . Details of the parameters associated with each identification
number are available as a dataframe by using the function **e2e_get_parmdoc()**.

\

The input distribution of parameter values to the Monte Carlo process is drawn from a random uniform distribution with limits specified in
the CredIntSim_SD control file for the model setup (located in a sub-folder of the /Parameters folder). This distribution is shown by
a red box and whisker at the bottom of the plot. Given the random uniform input we expect the quartiles (limits of the box) to be symetrical
and located mid-way between zero and the upper and lower extremes. Vertical red lines show the expected limits of the quartiles boxes if 
model results were completely insensitive to individual parameter values.

\

The extent to which individual parameter distributions deviate from the random uniform input is an indication of their sensitivity in the model. 
Parameters whose distributions are tightly confined around zero (the baseline value) are highly sensitive.

\

For some parameters, in particular the preference parameters, their credible distributions may span a wider range than the inputs. This may seem
unexpected, but arises because within each feeding guild the preference parameters are not independent of each other. The preferences within each guild
must sum to 1. Hence, during the Monte Carlo process, after drawing new values of the preference values they are all rescaled to sum to 1, which may mean that
some of them will have been varied by more than the original coefficient of variation of the input random uniform. 

\

For details of how the distribution of credible output values from StrathE2E are calculated see the help information for the e2e_run_mc() function.

\

```{r sensmcplot1, fig.width=7, fig.height=7, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package:
   model <- e2e_read("North_Sea", "1970-1999")
#
# Either run the functions e2e_run_sens() to generate some results, or use the internal example data
#
# Plot the ecology model parameter distributions:
   e2e_plot_sens_mc(model, selection="SENS", use.example=TRUE)  
```

**Figure 14**. Sensitivity analysis plot for the 1970-1999 North Sea model.

\

```{r sensmcplot2, fig.width=7, fig.height=9, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package:
   model <- e2e_read("North_Sea", "1970-1999")
#
# Either run the functions e2e_run_mc() to generate some results, or use the internal example data
#
# Plot the ecology model parameter distributions:
   e2e_plot_sens_mc(model, selection="MC", use.example=TRUE)  
```

**Figure 15**. Parameter likelihood distributions from a baseline mode Monte Carlo analysis with the 1970-1999 North Sea model.

\

\

# 7 Comparing model runs and observations

Comparing model results with observational data is a key part of the process of validation. Similarly, comparing outputs from different model runs is a 
key aspect of experimenting with the model.

\

The extensive model outputs provide scope for many creative ways of visualizing the relationships between models and observations and between models. Section 4 
provides guidance of accessing the model outputs for users to develop their own graphics. This section describes a set of standard functions for displaying model and observational 
data that have proved useful.

\

## 7.1 Comparing model results with obseravtions **e2e_compare_obs()**

Box and whisker plots comparing annual or monthly observational data with corresponding model outputs.

**Table 28**. Arguments of the function **e2e_compare_obs()**.

| Argument      | Description                                                                                                                    |
| ------------- | ------------------------------------------------------------------------------------------------------------------------------ |
| model         | R-list object generated by the e2e_read() function which defines the model configuration                                       |
| selection     | Text string from a list to select comparison with annual or monthly observations. Select from: "ANNUAL", "MONTHLY". Remember to include the phrase within "" quotes |
| ci.data       | Logical. If TRUE plot credible intervals around model results based on Monte Carlo simulation with the e2e_run_mc() function (default=FALSE)                      |
| use.saved     | Logical. If TRUE use data from a prior user-defined run held as csv files data in the current results folder (default=FALSE) |
| use.example   | Logical. If TRUE use pre-computed example data from the internal North Sea model rather than user-generated data (default=FALSE) |
| results       | R-list object of model output generated by the e2e_run() function (default=NULL) |


\

Generate a multi-panel set of box and whisker diagrams comparing annual or monthly averaged or integrated observational data on the state of the
ecosystem with corresponding model outputs. Each panel displays a different category of observational and model data. 

\

The observational data are read from the file 'annual_target_data*.csv' or 'monthly_target_data*.csv'in the /Target_data folder of the model setup defined by a e2e_read() function call.
Column 3 of annual_target_data_* (header "Use_1.yes_0.no") is a flag to set whether any given row is used in calculating the likelihood of the observed data given the model setup by functions 
such as e2e_optimize_eco(). Un-used rows of data are omitted from the box and whisker plotting panels.

\

Arguments determine the source of model data for comparison with the observations. These can be outputs from a single model run with data held in memory as a list object or in a saved csv file, or from 
a Monte Carlo simulation (using the function e2e_run_mc()) to estimate credible intervals of model outputs. Generation of credible interval data is a long computing task, so 
example data for the North Sea model provided with the package for illustration.

\

In each plot panel, the distribution of observation data is shown by a black box plot (box spans 50% of obervations, whiskers span 99%, median indicated by a black tick mark). Corresponding data 
derived from model outputs are always shown in red - either a red tick mark for data from a single model run, or a red boxplot for likelihood weighted data from a Monte Carlo analysis.

\

```{r compareobs1, fig.width=7, fig.height=7, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package, run, and generate a plot:
# Note that these are the observational data that were used as the target for optimizing the model parameters
model <- e2e_read("North_Sea", "1970-1999")
results <- e2e_run(model, nyears=3)
e2e_compare_obs(selection="ANNUAL", model, results=results)
```

**Figure 16**. Comparison of 1970-1999 annual integrated or averaged observational data from the North Sea (black box-plots) with corresponding values derived from a single run of the model (red symbols). 
For the observations, the box spans 50% of the distribution of observations and the whiskers 99%. The median is shown by a black tick mark.

\

```{r compareobs2, fig.width=7, fig.height=7, fig.align='center'}
# Compare results from the 1970-1999 version of the North Sea model with monthly averaged data:
# Note that these are the observational data that were NOT used for optimizing the model parameters
e2e_compare_obs(model, selection="MONTHLY",results=results)
```

**Figure 17**. Comparison of 1970-1999 monthly averaged observational data from the North Sea (black box-plots) with corresponding values derived from a single run of the model (red symbols). 
For the observations, the box spans 50% of the distribution of observations and the whiskers 99%. The median is shown by a black tick mark.

\

```{r compareobs3, fig.width=7, fig.height=7, fig.align='center'}
# As above, but plotting credible intervals of model data derived from a baseline mode Monte Carlo analysis:
model <- e2e_read("North_Sea", "1970-1999")
e2e_compare_obs(model, selection="ANNUAL",ci.data=TRUE,use.example=TRUE)
```

**Figure 18**. Comparison of 1970-1999 annual integrated or averaged observational data from the North Sea (black box-plots) with credible intervals of corresponding 
values derived from a Monte Carlo analysis of the model (red boxplots). 
For both observations and model outputs, the box spans 50% of the distribution of observations and the whiskers 99%. The median is shown by a tick mark.

\

```{r compareobs4, fig.width=7, fig.height=7, fig.align='center'}
e2e_compare_obs(selection="MONTHLY", model, ci.data=TRUE,use.example=TRUE)
```

**Figure 19**. Comparison of 1970-1999 monthly averaged observational data from the North Sea (black box-plots) with credible intervals of corresponding 
values derived from a Monte Carlo analysis of the model (red boxplots). 
For both observations and model outputs, the box spans 50% of the distribution of observations and the whiskers 99%. The median is shown by a tick mark.

\

\

## 7.2 Comparing baseline and scenario model results

Two functions are provided for comapring baseline and scenario model runs. The first relies on a boxplot format similar to that used for comapring model runs with 
observations. The second utilizes tornado barplots which are widely used in sensitivity analysis.

\

### 7.2.1 Comparing baseline and scenario model results using boxplots **e2e_compare_runs_box()**

Box and whisker plots comparing annual or monthly outputs from single or Monte Carlo baseline and scenario model runs.

**Table 29**. Arguments of the function **e2e_compare_runs_box()**.

| Argument      | Description                                                                                                                    |
| ------------- | ------------------------------------------------------------------------------------------------------------------------------ |
| selection     | Text string from a list to select comparison with annual or monthly observations. Select from: "ANNUAL", "MONTHLY". Remember to include the phrase within "" quotes |
| model1        | R-list object defining the baseline model configuration compiled by the e2e_read() function                                       |
| ci.data1      | Logical. If TRUE plot credible intervals around baseline model results based on Monte Carlo simulation with the e2e_run_mc() function (default=FALSE)                      |
| use.saved1    | Logical. If TRUE use baseline data from a prior user-defined run held as csv files data in the current results folder (default=FALSE) |
| use.example1  | Logical. If TRUE use pre-computed example data from the internal North Sea model as the baseline rather than user-generated data (default=FALSE) |
| results1      | R-list object of baseline model output generated by the e2e_run() function. Only needed if ci.data1=FALSE, use.saved1=FALSE and use.example1=FALSE. (Default=NULL) |
| model2        | R-list object defining the scenario model configuration compiled by the e2e_read() function                                       |
| ci.data2      | Logical. If TRUE plot credible intervals around scenario model results based on Monte Carlo simulation with the e2e_run_mc() function (default=FALSE)                      |
| use.saved2    | Logical. If TRUE use scenario data from a prior user-defined run held as csv files data in the current results folder (default=FALSE) |
| use.example2  | Logical. If TRUE use pre-computed example data from the internal North Sea model as the scenario rather than user-generated data (default=FALSE) |
| results2      | R-list object of scenario model output generated by the e2e_run() function. Only needed if ci.data1=FALSE, use.saved1=FALSE and use.example1=FALSE. (Default=NULL) |


Generate a multi-panel set of box and whisker diagrams comparing annual or monthly averaged or integrated model outputs
from single or Monte Carlo baseline (black) and scenario (red) model runs. Each panel of annual data displays a different category of model data; each panel of monthly data shows a different nutrient or plankton guild.

\

Arguments determine the source of model data for comparison. These can be outputs from 
a Monte Carlo simulations (using the function e2e_run_mc()) to estimate credible intervals of model outputs, or from single model runs using e2e_run(). Generation of credible interval data is a long computing task, so 
example data for the North Sea model provided with the package are available as an illustration.

\

In each plot panel, where credible intervals of model outputs (from Monte Carlo analysis) have been selected 
the box spans 50% of the lielihood distribution of data, whiskers span 99%, and the median is indicated by tick mark). Baseline model results are 
always shown in black, scenario results in red. 

\


```{r comparerunsbox1, fig.width=7, fig.height=7, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package and run for 3 years
m1 <- e2e_read("North_Sea", "1970-1999")
r1<-e2e_run(m1,nyears=3)
# Load the 2003-2013 version of the North Sea model supplied with the package and run for 3 years
m2 <- e2e_read("North_Sea", "2003-2013")
r2<-e2e_run(m2,nyears=3)
#Compare 1970-1999 as baseline (from example data with cred.int.), with 2003-2013 as sceanrio (from single model run)
e2e_compare_runs_box(selection="ANNUAL",  model1=m1, ci.data1=TRUE, use.example1=TRUE, model2=m2, ci.data2=FALSE, results2=r2)
```

**Figure 20**. Comparison of 1970-1999 annual integrated or averaged data from a single run of the North Sea (baseline model, black symbols) with corresponding 
values derived from a single run of the 2003-2013 model (scenario model, red symbols). 

\


```{r comparerunsbox2, fig.width=7, fig.height=7, fig.align='center'}
e2e_compare_runs_box(selection="MONTHLY", model1=m1, ci.data1=TRUE, use.example1=TRUE, model2=m2, ci.data2=FALSE, results2=r2)
```

**Figure 21**. Comparison of 1970-1999 monthly averaged data from a single run of the North Sea (baseline model, black symbols) with corresponding 
values derived from a single run of the 2003-2013 model (scenario model, red symbols). 

\

```{r comparerunsbox3, fig.width=7, fig.height=7, fig.align='center'}
#Compare 1970-1999 as baseline (from single model run), with 2003-2013 as sceanrio (from example data with cred.int.)
e2e_compare_runs_box(selection="ANNUAL",  model1=m1, ci.data1=FALSE, results1=r1, model2=m2, ci.data2=TRUE, use.example2=TRUE)
```

**Figure 22**. Comparison of 1970-1999 annual integrated or averaged data from a single run of the North Sea (baseline model, black symbols) with corresponding 
values derived from a Monte Carlo analysis of the 2003-2013 model (scenario model, red boxplots). The scenario model boxes span 50% of the distribution of 
likelihood weighted model results and the whiskers 99%. The median is shown by a tick mark.

\

```{r comparerunsbox4, fig.width=7, fig.height=7, fig.align='center'}
e2e_compare_runs_box(selection="MONTHLY", model1=m1, ci.data1=FALSE, results1=r1, model2=m2, ci.data2=TRUE, use.example2=TRUE)
```

**Figure 23**. Comparison of 1970-1999 monthly averaged data from a single run of the North Sea (baseline model, black symbols) with corresponding 
values derived from a Monte Carlo analysis of the 2003-2013 model (scenario model, red boxplots). The scenario model boxes span 50% of the distribution of 
likelihood weighted model results and the whiskers 99%. The median is shown by a tick mark.

\


```{r comparerunsbox5, fig.width=7, fig.height=7, fig.align='center'}
#Compare 1970-1999 as baseline (from example data with cred.int.), with 2003-2013 as sceanrio (from example data with cred.int.)
e2e_compare_runs_box(selection="ANNUAL", model1=m1, ci.data1=TRUE, use.example1=TRUE, model2=m2, ci.data2=TRUE, use.example2=TRUE )
```

**Figure 24**. Comparison of 1970-1999 annual integrated or averaged data from a Monte Carlo analysis of the North Sea (baseline model, black boxplots) with corresponding 
values derived from a Monte Carlo analysis of the 2003-2013 model (scenario model, red boxplots). In both cases the boxes span 50% of the distribution of 
likelihood weighted model results and the whiskers 99%. The median is shown by a tick mark.

\


```{r comparerunsbox6, fig.width=7, fig.height=7, fig.align='center'}
e2e_compare_runs_box(selection="MONTHLY", model1=m1, ci.data1=TRUE, use.example1=TRUE, model2=m2, ci.data2=TRUE, use.example2=TRUE )
```

**Figure 25**. Comparison of 1970-1999 monthly averaged data from a Monte Carlo analysis of the North Sea (baseline model, black boxplots) with corresponding 
values derived from a Monte Carlo analysis of the 2003-2013 model (scenario model, red boxplots). In both cases the boxes span 50% of the distribution of 
likelihood weighted model results and the whiskers 99%. The median is shown by a tick mark.

\

\

### 7.2.2 Comparing baseline and scenario model results using tornado barplots **e2e_compare_runs_bar()**

Create a tornado bar-plots of the differences between two model runs.

**Table 30**. Arguments of the function **e2e_compare_runs_bar()**.

| Argument      | Description                                                                                                                    |
| ------------- | ------------------------------------------------------------------------------------------------------------------------------ |
| selection     | Text string from a list identifying which type of data are to be plotted. Select from: "AAM", "CATCH", corresponding to annual average mass data, and catch data respectively (default = "AAM"). Remember to include the phrase within "" quotes |
| model1        | R-list object defining the baseline model configuration compiled by the e2e_read() function                                       |
| use.saved1    | Logical. If TRUE use baseline data from a prior user-defined run held as csv files data in the current results folder (default=FALSE) |
| results1      | R-list object of baseline model output generated by the e2e_run() function. Only needed if ci.data1=FALSE, use.saved1=FALSE and use.example1=FALSE. (Default=NULL) |
| model2        | R-list object defining the scenario model configuration compiled by the e2e_read() function                                       |
| use.saved2    | Logical. If TRUE use scenario data from a prior user-defined run held as csv files data in the current results folder (default=FALSE) |
| results2      | R-list object of scenario model output generated by the e2e_run() function. Only needed if ci.data1=FALSE, use.saved1=FALSE and use.example1=FALSE. (Default=NULL) |
| log.pc        | Value="LG" for data to be plotted on a log10 scale, value = "PC" for data to be plotted on a percentage difference scale (default = "PC") |
| zone          | Value = "O" for offshore, "I" for inshore, or "W" for whole model domain (all upper case) (default = "W") |
| bpmin         | Axis minimum for plot - i.e. the maximum NEGATIVE value of (scenario-baseline). Default = -50, given log.pc="PC" (percentage differences). Needs to be reset to e.g. -0.3 if log.pc="LG" (log scale) |
| bpmax         | Axis maximum for plot - i.e. the maximum POSITIVE value of (scenario-baseline). Default = +50, given log.pc="PC" (percentage differences). Needs to be reset to e.g. +0.3 if log.pc="LG" (log scale) |
| maintitle     | An optional descriptive text field (in quotes) to be added above the plot. Keep to 45 characters including spaces (default="") |

\

Create a tornado bar-plot diagram of the differences in either annual average masses of ecology model varibles, or annual fishery catches, between two
different runs of the StrathE2E model, referred to as baseline and sceanrio runs. 

\

A tornado plot is a horizontal barplot which shows the difference between baseline and sceanrio runs. Bars to the right indicate scenario values greater than the baseline, bars to the left indicate scenario values less than baseline.
In this function the difference between scenario and baseline data is calculated as either ((scenario - baseline)/baseline) and plotted on a linear (percentage) scale, or (scenario/baseline) and plotted on a log10 scale. In both cases, zero represents scenario = baseline.

\

The choice of whether to plot data on the masses of ecosystem components, or fishery catches, is determined by an argument setting in the function call. In either case two panels of tornado diagrams are plotted. In the case of mass data, the panels show
data on water column components of the ecoststem (upper panel) and seabed components of the system (lower panel). In the case of catch data, the upper panel shows landings, the lower shows discards. For mass data and landings, positive values (scenario > baseline) are
indicated by green bars to the right, negative values (scenario < baseline) by red bars to the left. For discards, positive values are in black, negative in grey. For any bars extending outsided the selected plotting range, the numeric value is displayed inside the relevant bar.

\

The function accepts data from either saved csv files created by the e2e_run() function, or directly from e2e_run() results objects.
The relevant csv files are: /results/Modelname/Variantname/ZONE_model_anav_biomass-*.csv where ZONE is INSHORE, OFFSHORE or WHOLEDOMAIN 
and * represents the model run identifier (model.ident) embedded in the R-list object created by the e2e_read() function.

\

In addition to generating graphics output the function returns a list object of the data presented in the plot. The object comprises two dataframes (changewater and changeseabed where annual average mass data are selected; changeland and changedisc where catch data are selected).
The first column in each dataframe is the proportional difference expressed on a log10 scale, the second column as a percentage.

\


```{r comparerunsbar1, fig.width=7, fig.height=7, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package :
base_model   <- e2e_read("North_Sea", "1970-1999",model.ident="baseline")
base_results <- e2e_run(base_model,nyears=3)
#
# Create a scenario run from 1970-1999 baseline:
scen1_model   <- base_model      # Copies the baseline configuration into a new model object
scen1_model$setup$model.ident <- "scenario1"
scen1_model$data$fleet.model$gear_mult[4] <- 0.5 # Gear 4 (Beam_Trawl_BT1+BT2) activity rate rescaled to 0.5*baseline
scen1_results <- e2e_run(scen1_model,nyears=30)
#
#Compare the annual average mass scenario data with the baseline
mdiff_results1 <- e2e_compare_runs_bar(selection="AAM",
                                       model1=NA, use.saved1=FALSE, results1=base_results,
				       model2=NA, use.saved2=FALSE, results2=scen1_results,
   	                               log.pc="PC", zone="W",
                 		       bpmin=(-40),bpmax=(+40),
       		                       maintitle="Beam Trawl activity reduced by half")
```

**Figure 26**. Tornado plot comparison of annual average masses of ecosystem components in a baseline model (North Sea, 1970-1999) with a scenario model in which beam trawl activity is reduced by half.
Red bars to the left indicate that the mass was lower in the scenario model than the baseline; green bars to the right indicate higher in the scenario model.

\

```{r comparerunsbar2, eval = TRUE}
str(mdiff_results1,max.level=1)   # View the structure of the returned list object
```


```{r comparerunsbar3, fig.width=7, fig.height=7, fig.align='center'}
#Create a second sceanario run
scen2_model   <- base_model     # Copies the baseline configuration into a new model object
scen2_model$setup$model.ident <- "scenario2"
scen2_model$data$fleet.model$gear_mult[1] <- 0.5 # Gear 1 (Pelagic_Trawl+Seine) activity rate rescaled to 0.5*baseline
scen2_results <- e2e_run(scen2_model,nyears=30)
#
#Compare the annual catches in the new scenario with the baseline previously saved in a csv files
mdiff_results2 <- e2e_compare_runs_bar(selection="CATCH",
			              model1=base_model, use.saved1=TRUE, results1=NA,
                     	              model2=NA,use.saved2=FALSE, results2=scen2_results,
                     	              log.pc="LG", zone="W",
                     	              bpmin=(-0.4),bpmax=(+0.6),
                     	              maintitle="Pelagic Trawl/Seine activity reduced by half")
```

**Figure 27**. Tornado plot comparison of annual catches and discards in a baseline model (North Sea, 1970-1999) with a scenario model in which pelagic trawl activity is reduced by half.
Upper panel (landings): red bars to the left indicate that the landed live weight was lower in the scenario model than the baseline; green bars to the right indicate higher in the scenario model.
Lower panel (discards): grey bars to the left indicate that the discarded live weight was lower in the scenario model than the baseline; black bars to the right indicate higher in the scenario model.

\

\

# 8 Generating fisheries yield curves

Fisheries yield curves are a key product that users will want to generate with the model. The yield curve relates stationary state catch of a guild (the catch that can be 
taken sustainable in the long term under stationary environmental and fishery conditions) to its harvest ratio. We expect the yield curve to be dome-shaped - at zero harvest 
ratio there is no catch; at very high harvest ratio the stock is depleted so there is little or no catch. At intermediate harvest ratios we expect a sustainable catch. 
The peak catch is referred to as the Maximum Sustainable Yield (MSY), and this is generated at a harvest ratio which is equivalent to the term $F_MSY$ in the fisheries literature. 
One of the most important tasks in fish stock assessment is to determine where the current fishing mortality rate of a species sits in relation to $F_MSY$.

\

Traditionally, MSY and $F_MSY$ are regarded as constant characteristics for each harvested species, or at least their potential for dependency on the status of other species is neglected. 
However, this is plainly not the case since changes in predation mortality must affect both MSY and $F_MSY$. This interdependency of yield properties of guilds is one the key features 
that StrathE2E was intended to resolve,

\

The function **e2e_run_ycurve()** runs a set of model scenarios defined by a vector of values of the harvest ratio multipler for 
either planktivorous or demersal fish. The baseline in each case is defined by a harvest ratio multipler vaklue of 1.0, i.e. the 'present-day' state of the fishery. For each muliplier 
in the vector, the model is run to a stationary state and the whole-domain annual average mass, annual landings and annual discards of planktivorous and demersal fish are extracted from the final year outputs. 
In addition the whole-domain annual average masses of all o fthe other living components of the food web are saved. The data are returned as a dataframe and saved to .csv files.

\

A separate function is provided to plot the yield curve data (**e2e_plot_ycurve())**.

\

## 8.1 Generating yield data **e2e_run_ycurve()**

Run a set of models to generate fishery yield curve data for either planktivorous or demersal fish.

**Table 31**. Arguments of the functions **e2e_run_ycurve()**

| Argument    | Description                                                                                                                    |
| ----------- | ------------------------------------------------------------------------------------------------------------------------------ |
| model       | R-list object defining a baseline model and configuration compiled by the e2e_read() function                                  |
| selection   | Text string from a list identifying the fish guild for which a yield curve is to be generated. Select from: "PLANKTIV", "DEMERSAL", Remember to include the phrase within "" quotes |
| nyears      | Number of years to run each model scenario (default = 50)                                                                      |
| HRvector    | A vector of ascending order unique multiplier values to be applied to the baseline fish harvest ratio of the guild to be analysed. The values do not have to be evenly spaced (default = c(0,0.5,1.0,1.5,2.0,2.5,3.0)) |
| HRfixed     | Single value of the multiplier to be applied as the alternative (planktivorous or demersal) fish harvest ratio for all runs (default=1.0) |
| csv.output  | Logical. If FALSE then disable writing of csv output files - useful for testing (default=TRUE) |

\

Perform a set of StrathE2E model runs along a sequence of values of either planktivorous or demersal fish harvest ratio multiplier, saving the annual average whole-domain biomass, annual landings and annual discards of all exploitable food web guilds from each run 
plus the annual average biomasses of all the other living components of the food web. 

\

The baseline for the sequence of runs (harvest ratio multiplier = 1.0) is a model name and variant as loaded by the e2e_read() function.

\

The planktivorous or demersal fish yield curve can be generated for a given fixed setting of the other (demersal or planktivorous) fish harvest ratio multiplier (default = 1.0). All other
conditions are held constant as in the baseline model configuration.

\

The yield curve represents the catch that would be generated from the stationary state of the model attained with long-term repeating annual cycles of all driving data. Hence it is important that each
simulation is run for long enough that the model attains its stationary state, which may be some distance from the baseline model initial conditions.
It is recommended that each run is at least 50 years.

\

The data on annual average biomass and annual integrated catches stored in the returned dataframe object (and optionally a csv file) can subsequently be plotted using the function e2e_plot_ycurve(). Users can easily plot any of the other 
saved data using their own plotting code.

\

If csv output is selected then the resulting files in the current user results folder have names Yield_curve_data_PFHRmult-*.csv or Yield_curve_data_DFHRmult-*.csv, depending on the 'selection' argument, and * represents the model.ident text set in the prior e2e_read() function call.

\

```{r runycurve, eval = FALSE}
# Load the 1970-1999 version of the North Sea model supplied with the package :
model <- e2e_read("North_Sea", "1970-1999",model.ident="70-99base")
#
# In this illustrative example the StrathE2E() model is run for only 5 years to enable quick return of results. In a real simulation nyear would be at least 50.
# This example illustrates that the vector of planktivorous fish harvest ratio multiplers does not have to be evenly spaced.
hr=c(0,0.5,0.75,1.0,1.25,2.0,3.0) 
pf_yield_data <- e2e_run_ycurve(model,selection="PLANKTIV", nyears=5,HRvector=hr,HRfixed=1,csv.output=TRUE)
#
# View the column names of the results dataframe:
names(pf_yield_data)
#
```

\

\

## 8.2 Plotting yield data **e2e_plot_ycurve()**

Plot fishery yield curve data for planktivorous or demersal fish.

**Table 32**. Arguments of the functions **e2e_plot_ycurve()**

| Argument     | Description                                                                                                                          |
| ------------ | ------------------------------------------------------------------------------------------------------------------------------------ |
| model        | R-list object defining the baseline model configuration used to generate the data and compiled by the e2e_read() function            |
| selection    | Text string from a list identifying the fish guild for which a yield curve is to be generated. Select from: "PLANKTIV", "DEMERSAL". Remember to include the phrase within "" quotes |
| use.saved    | Logical. If TRUE use data from a prior user-defined run held as csv files data in the current results folder (default=FALSE) |
| use.example  | Logical. If TRUE use pre-computed example data from the internal North Sea model rather than user-generated data (default=FALSE) |
| results      | Dataframe generated by the function e2e_run_ycurve(). Only needed if use.saved2=FALSE and use.example2=FALSE. (Default=NULL) |
| title        | Optional free text (enclosed in "") to be added as a header for the plot (default = "") |

\

Plot planktivorous or demersal fish yield curve data generated by the function e2e_run_ycurve(). 

\

In the function e2e_run_ycurve(), the baseline for the sequence of runs (harvest ratio multiplier = 1.0) is a model name and variant as loaded by the e2e_read() function.
The function then generates a set of biomass, landings and discards data for multiples of the target fish (planktivorous or demersal) harvest ratios relative to this baseline. This is done for a given value of
the alternative (demersal or planktivorous) harvest ratio (also a multiple of the the baseline). 

\

This function plots two graphs - the annual average fish biomass in the whole model domain (mMN.m^-2^) as a function of harvest ratio multiplier, and the yield curve, 
i.e. the annual catch (and discards) (mMN.m^-2^.y^-1^ from the whole model domain) as functions of the multiplier.

\

The yield curve represents the catch that would be generated from the stationary state of the model attained with long-term repeating annual cycles of all driving 
data including fishing.

\

Arguments for this function permit the input data to be drawn from an existing data object generated by the function e2e_run_ycurve(), a previously generated csv 
file, or example data provided with the package for versions of the internal North Sea models.

\

The function returns a dataframe containing the data which have been plotted plus the annual average mass density and where appropritae the landings and discards
of all the other living food web components in the model.

\


```{r plotycurve1, fig.width=7, fig.height=7, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package and generate a yield data object:
# WARNING - this example will take about several minutes to run
model <- e2e_read("North_Sea", "1970-1999")
pfhr=c(0,0.5,0.75,1.0,1.25,1.5,2.0,2.5,3.0) 
pf_yield_data <- e2e_run_ycurve(model,selection="PLANKTIV",nyears=50,HRvector=pfhr,HRfixed=1,csv.output=FALSE)
data <- e2e_plot_ycurve(model,selection="PLANKTIV", results=pf_yield_data,
                        title="Planktivorous yield with baseline demersal fishing")
```

\

```{r plotycurve2, fig.width=7, fig.height=7, fig.align='center'}
# Using examaple data generated with selection="PLANKTIV" ...
# Plot example data for one of the North Sea model versions internal to the package
model <- e2e_read("North_Sea", "1970-1999")
pf_yield_data<-e2e_plot_ycurve(model, selection="PLANKTIV", use.example=TRUE)
```
**Figure 28**. Yield data plot for planktivoirous fish in the 1970-1999 North Sea model with baseline demersal harvest ratios. Upper panel, annual average biomass of 
planktivorous fish (mMN.m^-2^ in the whole modle domain) as a function of harvest ratio. Lower panel: catch divided into landings and discards of planktivorous fish (mMN.m^-2^.y^-1^) as a function of harvest ratio. 

\

```{r plotycurve3, fig.width=7, fig.height=7, fig.align='center'}
# Users can then plot other biomass, landings and discards data in the results object by, for example:
par(mfrow=c(2,1))
par(mar=c(3.2,5,2,0.8))
ym<-1.1*max(pf_yield_data$Cetaceanbiom)
plot(pf_yield_data$PlankFishHRmult,pf_yield_data$Cetaceanbiom,ylim=c(0,ym),type="l",lwd=3,yaxt="n",xaxt="n",ann=FALSE)
   abline(v=1,lty="dashed")
   axis(side=1,las=1,cex.axis=0.9)
   axis(side=2,las=1,cex.axis=0.9)
   mtext("Planktiv. fish harvest ratio multiplier",cex=1,side=1,line=2)
   mtext("Cetacean biomass",cex=1,side=2,line=3.5)
   mtext(bquote("mMN.m"^-2),cex=0.7,side=3,line=-0.05,adj=-0.18)
ym<-1.1*max(pf_yield_data$Cetaceandisc)
plot(pf_yield_data$PlankFishHRmult,pf_yield_data$Cetaceandisc,ylim=c(0,ym),type="l",lwd=3,yaxt="n",xaxt="n",ann=FALSE)
   abline(v=1,lty="dashed")
   axis(side=1,las=1,cex.axis=0.9)
   axis(side=2,las=1,cex.axis=0.9)
   mtext("Planktiv. fish harvest ratio multiplier",cex=1,side=1,line=2)
   mtext("Cetacean by-catch",cex=1,side=2,line=3.5)
   mtext(bquote("mMN.m"^-2 ~ ".y"^-1),cex=0.7,side=3,line=-0.05,adj=-0.18)
```

**Figure 29**. Example of a possible user-generated yield data plot for cetaceans in the 1970-1999 North Sea model as a function of planktivorous harvest ratio, with baseline demersal harvest ratios. Upper panel, annual average biomass of 
cetaceans (mMN.m^-2^ in the whole modle domain) as a function of planktivorous fish harvest ratio. Lower panel: catch divided into targeted landings and by-catch (discards) of cetaceans (mMN.m^-2^.y^-1^) as a function of harvest ratio. 

\

```{r plotycurve4, fig.width=7, fig.height=7, fig.align='center'}
# Using example data generated with selection="DEMERSAL"...
model <- e2e_read("North_Sea", "1970-1999")
df_yield_data<-e2e_plot_ycurve(model, selection="DEMERSAL", use.example=TRUE)
```

**Figure 30**. Yield data plot for demersal fish in the 1970-1999 North Sea model with baseline planktivorous harvest ratios. Upper panel, annual average biomass of 
demersal fish (mMN.m^-2^ in the whole modle domain) as a function of harvest ratio. Lower panel: catch divided into landings and discards of demersal fish (mMN.m^-2^.y^-1^) as a function of harvest ratio. 

\

```{r plotycurve5, fig.width=7, fig.height=7, fig.align='center'}
# Users can then plot other biomass, landings and discards data in the results object by, for example:
par(mfrow=c(2,1))
par(mar=c(3.2,5,2,0.8))
ym<-1.1*max(df_yield_data$Pinnipedbiom)
plot(df_yield_data$DemFishHRmult,df_yield_data$Pinnipedbiom,ylim=c(0,ym),type="l",lwd=3,yaxt="n",xaxt="n",ann=FALSE)
    abline(v=1,lty="dashed")
    axis(side=1,las=1,cex.axis=0.9)
    axis(side=2,las=1,cex.axis=0.9)
    mtext("Demersal fish harvest ratio multiplier",cex=1,side=1,line=2)
    mtext("Pinniped biomass",cex=1,side=2,line=3.7)
    mtext(bquote("mMN.m"^-2),cex=0.7,side=3,line=-0.05,adj=-0.18)
ym<-1.1*max(df_yield_data$Pinnipeddisc)
plot(df_yield_data$DemFishHRmult,df_yield_data$Pinnipeddisc,ylim=c(0,ym),type="l",lwd=3,yaxt="n",xaxt="n",ann=FALSE)
    abline(v=1,lty="dashed")
    axis(side=1,las=1,cex.axis=0.9)
    axis(side=2,las=1,cex.axis=0.9)
    mtext("Demersal fish harvest ratio multiplier",cex=1,side=1,line=2)
    mtext("Pinniped by-catch",cex=1,side=2,line=3.7)
    mtext(bquote("mMN.m"^-2 ~ ".y"^-1),cex=0.7,side=3,line=-0.05,adj=-0.18)
```

**Figure 30**. Example of a possible user-generated yield data plot for pinnipeds in the 1970-1999 North Sea model as a function of demersal harvest ratio, with baseline planktivorous harvest ratios. Upper panel, annual average biomass of 
pinnipeds (mMN.m^-2^ in the whole modle domain) as a function of demersal fish harvest ratio. Lower panel: by-catch (discards) of pinnipeds (mMN.m^-2^.y^-1^) as a function of harvest ratio. In this case there are no targeted landings.

\

\

# 9 Plotting functions for visualizing model inputs

There are two main functions for plotting the input data to the model - **e2e_plot_edrivers()** and e2e_plot_fdrivers()** - these 
visualize the environmental inputs and fishing-related inputs respectively.

\

## 9.1 Visualizing environmental drivers **e2e_plot_edrivers()**

Plot climatological year of environmental driving data.

**Table 33**. Arguments of the functions **e2e_plot_edrivers()**

| Argument     | Description                                                                                     |
| ------------ | ----------------------------------------------------------------------------------------------- |
| model        | R-list object defining a model configuration and compiled by the e2e_read() function            |

\

Create a multi-panel page of time series plots of climatological annual cycles of driving data as provided in the input csv files.

\

The variables to be plotted and their units are shown in Table 34.

\

**Table 34**. Environmental input data dispalayed in each panel of the plot produced by **e2e_plot_edrivers()**


| Data                                          | Units                 | Comments                                                         |
| --------------------------------------------- | --------------------- | ---------------------------------------------------------------- |
| Sea surface irradiance                        | $\mu$E.m^-2^.d^-1^            |                                                                  |
| Suspended particulate matter                  | g.m^-3^                 |                                                                  |
| Temperature                                   | deg-C                 |
| Vertical diffusivity gradient                 | m.d^-1^                 | derived from the vertical diffusivity (m^2^.s^-1^) and mixing length scale (m) |
| External inflows                              | m^3^ per m^2^ sea surface of model domain |derived from proportion input per layer volume, layer thicknesses and areas |
| River discharge                               | m^3^ per m^2^ sea surface of model domain  | derived from proportion input to inshore volume, and inshore layer thickness and area |
| Inshore significant wave height               | m                     |                                                                  |
| Proportion of seabed disturbed                | d^-1^                   | aggregated over the three sediment classes in each zone          |
| External boundary nitrate concentration       | mMN.m^-3^               |                                                                |
| External boundary ammonia concentration       | mMN.m^-3^               |                                                                |
| External boundary phytoplankton concentration | mMN.m^-3^               |                                                                |
| External boundary detritus concentration      | mMN.m^-3^               |                                                                |
| River nitrate concentration                   | mMN.m^-3^               |                                                                |
| River ammonia concentration                   | mMN.m^-3^               |                                                                |
| Atmospheric nitrate deposition flux           | mMN.m^-2^.d^-1^           | area averaged flux density over each zone                       |
| Atmospheric ammonia deposition flux           | mMN.m^-2^.d^-1^           | area averaged flux density over each zone                       |

\

```{r plotenvdrivers, fig.width=7, fig.height=7, fig.align='center'}
model<-e2e_read(model.name="North_Sea", model.variant="2003-2013")
e2e_plot_edrivers(model)
```

**Figure 32**. Plot of climatological annual cycles of the environmental drivers of the 2003-2013 North Sea model


\

\

## 9.2 Visualizing fishing-related drivers **e2e_plot_fdrivers()**

Plot shaded-grid-cell diagrams of the distributions of fishing-related drivers across gears, guilds and seabed habitats.

**Table 35**. Arguments of the functions **e2e_plot_fdrivers()**

| Argument     | Description                                                                                     |
| ------------ | ----------------------------------------------------------------------------------------------- |
| model        | R-list object defining a model configuration and compiled by the e2e_read() function            |
| selection    | Text string from a list identifying the class of fishing-related driving data to be plotted. Select from: "ACTIVITY", "ABRASION", "HARVESTR", "DISCARDS", "OFFAL". Remember to include the phrase within "" quotes |

\

Select from a list of fishing-related drivers of the model (gear activity rates, seabed abrasion rates, harvest ratios, discard rates and offal generation rates), and 
generate a shaded-grid-cell diagram showing the distribution across gears, guilds and seabed habitats.

\

In each case, regardless of the vaviable selected, the function returns both a graphical image and a list object containing the matix of data plotted
and vectors of the axis labels, which uses can use to generate their own plots if desired.

\

**Details relating to fishing activity distributions (selection = "ACTIVITY"):**

Plot a matrix of seabed habitats vs fishing gears with each cell shaded to indicate activity density (white = 0, purple = high).

\

The spatial distribution of fishing gear activity in the model is defined by two input data sets:
* Vector of whole domain activity density of each fishing gear (s.d^-1^ per m^2^ of whole model domain)
* Matrix of the proportional distribution of whole domain activity density of each gear (rows) across seabed habitats (columns)

\

The activity density of each gear in each habitat is then obtained by multiplying the whole domain activity density into the proportional
distribution matrix, and dividing by the area-proportions of habitats in the domain. The units of habitat-specific activity density
are then (s.d^-1^.m^-2^).

\

The vector of area-proportions of each habitat in the model domain is part of the moldel configuration parameter set.

\

The calculation take account of any activity multipliers set in the csv inputs to be applied to the activity density.

\

When interpreting patterns of activity density, bear in mind that the impact (harvest ratio, seabed abrasion) per unit activity can vary considerably between gears.

\

```{r plotfdrivers1, fig.width=7, fig.height=6, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package, run, and generate a plot:
model <- e2e_read("North_Sea", "2003-2013")
plotted_data<-e2e_plot_fdrivers(model, selection="ACTIVITY")
```

**Figure 33**. Distribution of gear activity rates across habitats in the 2003-2013 North Sea model.

\

**Details relating to seabed abrasion rate distributions (selection = "ABRASION"):**


The function plots a matrix of seabed habitats vs fishing gears with each cell shaded to indicate abrasion area ratio (log_10_ scaled proportion of habitat area abraded per day; white = 0, purple = high).

\

The spatial distribution of seabed abrasion in the model is defined by three input data sets:
* Vector of whole domain activity density of each fishing gear (s.d^-1^ per m^2^ of whole model domain)
* Matrix of the proportional distribution of whole domain activity density of each gear (rows) across seabed habitats (columns)
* Vector of seabed abrasion rate of each fishing gear (m^2^.s^-1^)

\

The activity density of each gear in each habitat is then obtained by multiplying the whole domain activity density into the proportional
distribution matrix, and dividing by the area-proportions of habitats in the domain. The units of habitat-specific activity density
are then (s.d^-1^.m^-2^).

\

The product of activity density per gear and habitat and gear abrasion rate is then the abrasion area ratio, or abrasion rate (m^2^.m^-2^.d^-1^)

\

The vector of area-proportions of each habitat in the model domain is part of the moldel configuration parameter set.

\

The calculation take account of any activity multipliers set in the csv inputs to be applied to the activity density.

\

```{r plotfdrivers2, fig.width=7, fig.height=6, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package, run, and generate a plot:
model <- e2e_read("North_Sea", "2003-2013")
plotted_data<-e2e_plot_fdrivers(model, selection="ABRASION")
```

**Figure 34**. Distribution of abrasion area ratio across habitats in the 2003-2013 North Sea model.

\

**Details relating to harvest ratio distributions (selection = "HARVSTR"):**

The function plots a matrix of the values of harvest ratio (proportion of exploitable biomass of each guild capatured per day) for each guild in the ecology model (columns) arising from the activity of each gear (rows). 
Hence each row represents the selectivity pattern of each gear with respect to guilds.
Cells in the matrix are shaded to indicate log_10_ transformed harvest ratio (colour gradient: white = 0, purple = high).

\

Harvest ratio is calculated by the fishing fleet model and piped into the ecology model. The fleet model takes inputs of activity density by each gear,
and a matrix of fishing power, to calculate effort (effort = activity * power). Harvest ratio is assumed to be lineally related to effort by a scaling
coefficient which is input as a parameter to the fleet model.

\

Note that the visualization of harvest ratio generated by this function is based on the csv input data to the fleet model (including any multipliers to be applied to
either fishing activity or harvest ratio scaling parameters). However, it does not reflect any effects on harvest ratio arising from discarding scenarios since
these are dynamic and generated at run-time within the ecology model.

\

```{r plotfdrivers3, fig.width=7, fig.height=6, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package, run, and generate a plot:
model <- e2e_read("North_Sea", "2003-2013")
plotted_data<-e2e_plot_fdrivers(model, selection="HARVESTR")
```

**Figure 35**. Distribution of harvest ratios by gears across guilds in the 2003-2013 North Sea model.

\

**Details relating to discard rate distributions (selection = "DISCARDS"):**

The function plots a matrix of the values of discard rate (proportion by weight of catch discarded at sea prior to any at-sea processing (evisceration)) for each guild in the ecology model (columns) arising from the activity of each gear (rows). 
Cells in the matrix are shaded to indicate discard rate (range 0 - max; colour gradient: white = 0, purple = high) on a linear scale.

\

It is possible (and allowable) for the discard rate of a gear-guild combination to be set as a positive number in the 'fishing_discard_parameters' csv input file, but nevertheless the
the catching power settings in the corresponding 'fishing_power_parameters' file to be set to zero - in other words the given gear does not actually catch the given guild. In such cases this
function resets the discard rate to NA for plotting purposes.

\

Note that the visualization of discard rate generated by this function is based on the csv input data to the fleet model. 
Hence, it does not reflect any effects on discard rate arising from discarding scenarios configured in the fleet model input 
since these are dynamic and generated during run-time within the ecology model.

\

```{r plotfdrivers4 , fig.width=7, fig.height=6, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package, run, and generate a plot:
model <- e2e_read("North_Sea", "2003-2013")
plotted_data<-e2e_plot_fdrivers(model, selection="DISCARDS")
```

**Figure 36**. Distribution of discard rates by gears across guilds in the 2003-2013 North Sea model.

\

**Details relating to offal generation rate distributions (selection = "OFFAL"):**

The function plots a matrix of the values of offal generation rate (i.e. the proportion by weight of total catch which is returned to the sea as viscera after processing) for each guild (columns) by each gear (rows). 
Cells in the matrix are shaded to indicate offal generation rate (range 0 - max; colour gradient: white = 0, purple = high) on a linear scale.

\

The offal generation rate is a quantity derived from three input variables to the fishing fleet model:
* Matrix of the discard rate (proportion of catch weight discarded without being processed) for each guild and gear; (D); (input file '/Parameters/fishing_discard_parameters*.csv')
* Matrix of the processing-at-sea rate (proportion of retained catch which is processed) for each guild and gear; (P); (input file '/Parameters/fishing_processing_at_sea_parameters*.csv')
* Proportion by weight of viscera in the retained catch (assumed constant across all guilds); (V); (value in input file '/Parameters/fishing_fleet_parameters*.csv')

The offal generation rate is then given by (1 - D).P.V

\

It is possible (and allowable) for the discard rate and processing-at-sea rate of a gear-guild combination to be set as positive numbers in the relevant input files, but nevertheless the
the catching power settings in the corresponding 'fishing_power_parameters' file to be set to zero - in other words the given gear does not actually catch the given guild. In such cases this
function resets the discard rate and processing-at-sea rate to NA for plotting purposes.

\

Note that the visualization of offal rate generated by this function is based on the csv input data to the fleet model. 
Hence, it does not reflect any effects on discard rate arising from discarding scenarios configured in the fleet model input 
since these are dynamic and generated during run-time within the ecology model.

\

```{r plotfdrivers5 , fig.width=7, fig.height=6, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package, run, and generate a plot:
model <- e2e_read("North_Sea", "2003-2013")
plotted_data<-e2e_plot_fdrivers(model, selection="OFFAL")
```

**Figure 37**. Distribution of offal generation rates by gears across guilds in the 2003-2013 North Sea model.

\

\

# 10 Plotting functions for visualizing model outputs from the final year of a run

The model generates a wide range of outputs from each run. Full details are given in the Technical Manual, and guidance on 
how to unpack the list object returned at the end of a run is given in Section 4. There is considerable scope for creativity
in generating graphics to visualize the data.

\

The package contains 5 main functions for visualizing model output. They do not provide comprehensive coverage of all the available 
outputs, more an entry point for users.

\

Variables to be plotted are selected by argument settings in the functions. In addition, the functions are configured to plot either 
data from single-runs of the model, or outputs from Monte Carlo simulations showing credible intervals of model outputs.

\

All of the functions in this section (**Table 36**) show outputs from the final year of any model runs, not transients. The philosophy of the model 
is focussed on stationary state solutions, and less on transient behaviours. This is because transients are dependent on initial conditions which are 
unknown in most cases, whereas stationary states are independent of initial conditions.

\

**Table 36**. Functions available for visualizing model outputs from the final year of a run

| Function             | Description                                                    |
| -------------------- | -------------------------------------------------------------- |
| e2e_plot_eco()       | Plot annual cycles of ecology model variables                  |
| e2e_plot_migration() | Plot annual cycles of active migration fluxes                  |
| e2e_plot_catch()     | Plot distributions of catches by gear and guild                |
| e2e_plot_trophic()   | Plot mean trophic level and omnivory indices                   |
| e2e_plot_biomass()   | Plot zonal distributions of annual average biomass densities   |

\

## 10.1 Plotting annual cycles of ecology model variables **e2e_plot_eco()**

Plot daily data on ecological outputs from the model over the final year of a run, optionally with credible intervals.

**Table 37**. Arguments of the functions **e2e_plot_eco()**

| Argument      | Description                                                                                                                    |
| ------------- | ------------------------------------------------------------------------------------------------------------------------------ |
| model         | R-list object defining the model configuration compiled by the e2e_read() function                                             |
| selection     | Text string from a list identifying the group of model output variables to be plotted. Select from: "NUT_PHYT", "SEDIMENT", "ZOOPLANKTON", "FISH", "BENTHOS", "PREDATORS", "CORP_DISC", "MACROPHYTE", Remember to include the phrase within "" quotes |
| ci.data       | Logical. If TRUE plot credible intervals around baseline model results based on Monte Carlo simulation with the e2e_run_mc() function (default=FALSE)                      |
| use.saved     | Logical. If TRUE use baseline data from a prior user-defined run held as csv files data in the current results folder (default=FALSE) |
| use.example   | Logical. If TRUE use pre-computed example data from the internal North Sea model as the baseline rather than user-generated data (default=FALSE) |
| results       | R-list object of model output generated by the e2e_run() function. Only needed if ci.data=FALSE, use.saved=FALSE and use.example=FALSE. (Default=NULL)

\

Generate a multi-panel set of one-year time series plots of selected outputs on the concentrations of ecology model state variables in each spatial zone (**Table 38**). The default is to plot data from a single model run but, if available, 
credible intervals of model output from a Monte Carlo analysis can be plotted instead.

\

**Table 38**. Variables plotted (and units) given values of the argument "selection"

| selection   | Variables plotted                                             |    Units               |
| ----------- | ------------------------------------------------------------- | ---------------------- |
| NUT_PHYT    | Water column nitrate, ammonia, detritus and phytoplankton     | mMN.m^-3^              |
| SEDIMENT    | Sediment porewater nitrate and ammonia, detritus and corpses  | mMN.m^-3^ or gN.gDW^-1^  |
| ZOOPLANKTON | Omnivorous and carnivorous zooplankton                        | mMN.m^-2^              |
| FISH        | Planktivorous and demersal fish and fish larvae               | mMN.m^-2^              |
| BENTHOS     | Susp/dep and carn/scav benthos and benthos larvae             | mMN.m^-2^              |
| PREDATORS   | Birds, pinnipeds, cetaceans and migratory fish                | mMN.m^-2^              |
| CORP_DISC   | Corpses and discards                                          | mMN.m^-2^              |
| MACROPHYTE  | Macrphytes and macrophyte debris                              | mMN.m^-2^              |

\

Arguments determine the source of model data to be plotted. These can be outputs from a single model run with data held in memory as a list object or in a saved csv file, or from 
a Monte Carlo simulation (using the function **e2e_run_mc()**) to estimate credible intervals of model outputs. Generation of credible interval data is a long computing task, so 
example data for the North Sea model provided with the package are available as an illustration.

\

If plotting of credible intervals is selected, results from the maximum likelihood model are shown by a red line. The median of the credible values distribution
is shown my a solid black line. A grey-shaded area indicates the 50% credible interval (spanning quartiles of the cumulative likelihood
of simulated values). Black dashed lines span the 99% credible interval.

\

To direct the graph output to a file rather than the screen, wrap the e2e_plot_eco() function call in a graphical device call:
Since the plot pages contain different numbers of panels the recommended width:height ratios are as follows:

\

**Table 39**. Suggested width:height ratios for graphical output to a file depending on the argument "selection"

| selection   | width | height|
| ----------- | ----- | ----- |
| NUT_PHYT    |  1.5  |   1   |
| SEDIMENT    | 0.67  |   1   |
| ZOOPLANKTON |    1  |   1   |
| FISH        |    2  |   1   |
| BENTHOS     |    2  |   1   |
| PREDATORS   |    2  |   1   |
| CORP_DISC   |    1  |   1   |
| MACROPHYTE  |    2  |   1   |

\

```{r ploteco1 , fig.width=7, fig.height=6, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package, run, and generate a plot:
model <- e2e_read("North_Sea", "1970-1999")
results <- e2e_run(model, nyears=3)
e2e_plot_eco(model, selection="NUT_PHYT",results=results)
```

**Figure 38**. Annual cycles of nutrients and phytoplankton in the final year of a 1970-1999 North Sea model.

\

```{r ploteco2 , fig.width=7, fig.height=7, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package and plot example credible interval data:
model <- e2e_read("North_Sea", "1970-1999")
e2e_plot_eco(model, selection="NUT_PHYT", ci.data=TRUE, use.example=TRUE)
```

**Figure 39**. Annual cycles of zooplankton with credible intervals in the final year of a 1970-1999 North Sea model.

\

Example of directing graphical output to a pdf file
```{r ploteco0 , eval=FALSE}
pdf("C:/Users/username/Documents/Foldername/plot.pdf",width=7,height=4.667)      # or jpeg("plot.jpg"), png("plot.png")
     e2e_plot_eco(model, "NUT_PHYT", ci.data=TRUE, use.example=TRUE)
dev.off()
```

\

\

## 10.2 Plotting annual cycles of active migration fluxes **e2e_plot_migration()**

Plot daily data on migration fluxes by actively mobile guilds during the final year of a run, optionally with credible intervals.

**Table 40**. Arguments of the functions **e2e_plot_migration()**

| Argument      | Description                                                                                                                    |
| ------------- | ------------------------------------------------------------------------------------------------------------------------------ |
| model         | R-list object defining the model configuration compiled by the e2e_read() function                                             |
| selection     | Text string from a list identifying the group of model output variables to be plotted. Select from: "NUT_PHYT", "SEDIMENT", "ZOOPLANKTON", "FISH", "BENTHOS", "PREDATORS", "CORP_DISC", "MACROPHYTE", Remember to include the phrase within "" quotes |
| ci.data       | Logical. If TRUE plot credible intervals around baseline model results based on Monte Carlo simulation with the e2e_run_mc() function (default=FALSE)                      |
| use.saved     | Logical. If TRUE use baseline data from a prior user-defined run held as csv files data in the current results folder (default=FALSE) |
| use.example   | Logical. If TRUE use pre-computed example data from the internal North Sea model as the baseline rather than user-generated data (default=FALSE) |
| results       | R-list object of model output generated by the e2e_run() function. Only needed if ci.data=FALSE, use.saved=FALSE and use.example=FALSE. (Default=NULL)

\

Generate a multi-panel set of one-year time series plots of the mass fluxes between inshore and offshore zones due to migration by actively mobile guilds 
in the ecology model: all three fish guilds, birds, pinnipeds and cetaceans. The default is to plot data from a single model run but if available, 
credible intervals of model output from a Monte Carlo analysis can be plotted instead.

\

Daily interval post-processed data from the Monte Carlo **e2e_run_mc()** function are stored in the file
/results/Modelname/Variantname/CredInt/CredInt_processed_daily_migrations-*.csv, where * represents the model run identifier
(model.ident) text embedded in the R-list object created by the **e2e_read()** function.

\

Arguments determine the source of model data to be plotted. These can be outputs from a single model run with data held in memory as a list object or 
in a saved csv file, or from a Monte Carlo simulation (using the function **e2e_run_mc()**) to estimate credible intervals of model outputs. Generation 
of credible interval data is a long computing task, so 
data for the North Sea model provided with the package are available as example data sets.

\

Each panel of the plot shows a time-series of the net flux densities (mMN.d^-1^ in the model domain as a whole, assumed to be 1 m^2^ sea surface area) of one
of the migratory guilds in the model (all three guilds of fish, birds, pinnipeds and cetaceans) between the inshore and offshore zones
of the model, over the final year of a run. These migration fluxes are the dynamic product of gradients in the ratio of food
concentration to predator concentration across the inshore-offshore boundary. Positive values of the net migration flux indicate net movement from the 
offshore to inshore zone. Negative values indicate net movement from inshore to offshore.

\

If plotting of credible intervals is selected, results from the maximum likelihood model are shown by a red line. The median of the credible values distribution
is shown my a solid black line. A grey-shaded area indicates the 50% credible interval (spanning quartiles of the cumulative likelihood
of simulated values). Black dashed lines span the 99% credible interval.

\

For details of how the distribution of credible output values from StrathE2E are calculated see help(e2e_run_mc).

```{r plotmig1 , fig.width=7, fig.height=7, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package, run, and generate a plot:
model <- e2e_read("North_Sea", "1970-1999")
results <- e2e_run(model, nyears=3)
e2e_plot_migration(model, results=results)
```

**Figure 40**. Annual cycles of migration fluxes between the offshore and inshore zones in the final year of a 1970-1999 North Sea model.

\

```{r plotmig2 , fig.width=7, fig.height=7, fig.align='center'}
# For the same model as supplied with the package, plot the example data with credible intervals:
model <- e2e_read("North_Sea", "1970-1999")
e2e_plot_migration(model, ci.data=TRUE, use.example=TRUE)
```

**Figure 41**. Annual cycles of migration fluxes with credible intervals between the offshore and inshore zones in the final year of a 1970-1999 North Sea model.

\

\

## 10.3 Plotting distributions of catches by gear and guild **e2e_plot_catch()**

Plot the distribution and composition of annual catches across guild, zones and gears in the final year of a model run.

**Table 41**. Arguments of the functions **e2e_plot_catch()**

| Argument      | Description                                                                                                                    |
| ------------- | ------------------------------------------------------------------------------------------------------------------------------ |
| model         | R-list object defining the model configuration compiled by the e2e_read() function                                             |
| results       | List object of single-run model output generated by running the function e2e_run() function
| selection     | Text string from a list identifying the group of model output variables to be plotted. Select from: "BY_GUILD", "BY_GEAR". With the former, each panel represents a different guild; with the latter, each panel represnts a different gear. Remember to include the phrase within "" quotes |

\

Create stacked barplots of the distribution of offshore and inshore landings and discards across gears by guild, or across guilds by gears, in the final year of a model run.

\

Data in zonal landings and discards by guild and gear in the final year of a model run are generated as a standard output from the model, and saved both as csv files and in the results object returned by the e2e_run() function.
This function organises these data in two different ways to display as barplots.

\

The first display is a multi-panel plot in which each panel represents a different guild, and the bars show the zonal landings by each gear. The alternative display
has each panel as a different gear, and the bars show the landings and discards of each guild.

\

The unit of the displayed data are mMN.y^-1^ from the model domain as a whole, which is taken as being 1 m^2^.

\

```{r plotcatch1 , fig.width=7, fig.height=7, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package, run, and generate a plot:
model <- e2e_read("North_Sea", "1970-1999")
results <- e2e_run(model, nyears=3)
e2e_plot_catch(model, results, selection="BY_GEAR")
```

**Figure 42**. Annual landings and discards of each guild, by gear, in the final year of a run of the 1970-1999 North Sea model.

\

```{r plotcatch2 , fig.width=7, fig.height=7, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package, run, and generate a plot:
e2e_plot_catch(model, results, selection="BY_GUILD")
```

**Figure 43**. Annual landings and discards by each gear, by guild, in the final year of a run of the 1970-1999 North Sea model.

\

\

## 10.4 Plotting mean trophic level and omnivory indices **e2e_plot_trophic()**

Plot showing the annual mean trophic level index, and the omnivory index of each guild during the final year of a model run, optionally with credible intervals.

**Table 42**. Arguments of the functions **e2e_plot_trophic()**

| Argument      | Description                                                                                                                    |
| ------------- | ------------------------------------------------------------------------------------------------------------------------------ |
| model         | R-list object defining the model configuration compiled by the e2e_read() function                                             |
| ci.data       | Logical. If TRUE plot credible intervals around baseline model results based on Monte Carlo simulation with the e2e_run_mc() function (default=FALSE)                      |
| use.saved     | Logical. If TRUE use baseline data from a prior user-defined run held as csv files data in the current results folder (default=FALSE) |
| use.example   | Logical. If TRUE use pre-computed example data from the internal North Sea model as the baseline rather than user-generated data (default=FALSE) |
| results       | R-list object of model output generated by the e2e_run() function. Only needed if ci.data=FALSE, use.saved=FALSE and use.example=FALSE. (Default=NULL)

\

Generate a two-panel plot showing: (upper panel) the mean trophic level of each guild in the ecology model, and (lower panel) the omnivory index of each guild. The data are 
generated by the NetIndices package from a flow matrix of nutrient fluxes through, into and out of the ecosystem during the final year of a run. The data are generated
automatically as part of the output from every call of the **e2e_run()** function. The default is to plot data from a single model run but if available, 
credible intervals of model output from a Monte Carlo analysis can be plotted instead.

\

If credible intervals are plotted these are displayed as box-and-whiskers. The box spans 50% of the likelihood distribution of values and the whiskers 99%. The median is shown 
by a black tick mark and the maximum likelihhod model by a red tick mark.

\

Note the the NetINdices package assigns trophic level 1 to inorganic nutrients and detritus, to autotrophs and detritivores are ranked as trophic level 2.

Arguments determine the source of model data to be plotted. These can be outputs from a single model run with data held in memory as a list object or in a saved csv file, or from 
a Monte Carlo simulation (using the function **e2e_run_mc()**) to estimate credible intervals of model outputs. Generation of credible interval data is a long computing task, so 
example data for the North Sea model provided with the package are available as illustration.

\


```{r plottrophic1 , fig.width=7, fig.height=7, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package, run, and generate a plot:
model <- e2e_read("North_Sea", "1970-1999")
results <- e2e_run(model, nyears=3)
e2e_plot_trophic(model, results=results)
```

**Figure 44**. Annual mean trophic level and omnivory indices for the final year of a single run of the 1970-1999 North Sea model.

\

```{r plottrophic2 , fig.width=7, fig.height=7, fig.align='center'}
# For the same model as supplied with the package, plot the example data with credible intervals:
e2e_plot_trophic(model, ci.data=TRUE, use.example=TRUE)
```

**Figure 45**. Annual mean trophic level and omnivory indices with credible intervals, for the final year of the 1970-1999 North Sea model.

\

\

## 10.5 Plotting zonal distributions of annual average biomass densities **e2e_plot_biomass()**

Plot showing the annual annual average biomass densities of each guild in the inshore and offshore zones, optionally with credible intervals.

**Table 43**. Arguments of the functions **e2e_plot_biomass()**

| Argument      | Description                                                                                                                    |
| ------------- | ------------------------------------------------------------------------------------------------------------------------------ |
| model         | R-list object defining the model configuration compiled by the e2e_read() function                                             |
| ci.data       | Logical. If TRUE plot credible intervals around baseline model results based on Monte Carlo simulation with the e2e_run_mc() function (default=FALSE)                      |
| use.saved     | Logical. If TRUE use baseline data from a prior user-defined run held as csv files data in the current results folder (default=FALSE) |
| use.example   | Logical. If TRUE use pre-computed example data from the internal North Sea model as the baseline rather than user-generated data (default=FALSE) |
| results       | R-list object of model output generated by the e2e_run() function. Only needed if ci.data=FALSE, use.saved=FALSE and use.example=FALSE. (Default=NULL)

\

Generate plots showing the annual average biomass densities of each guild in the ecology model in the inshore (blue) and offshore (red) zones during the final year of a run. The default is to plot data from a single model run but if available, 
credible intervals of model output from a Monte Carlo analysis can also be plotted.

\

If credible intervals are plotted these are displayed as box-and-whiskers. The box spans 50% of the likelihood distribution of values and the whiskers 99%. The median is shown 
by a tick mark.

\

Note that in this plot the biomass are expressed as mMN.m^-2^, meaning that the mass in each zone has been scaled to the zonal sea surface area. So the data in
this plot are area densities and not mass.

\

Arguments determine the source of model data to be plotted. These can be outputs from a single model run with data held in memory as a list object or in a saved csv file, or from 
a Monte Carlo simulation (using the function **e2e_run_mc()**) to estimate credible intervals of model outputs. Generation of credible interval data is a long computing task, so 
example data for the North Sea model provided with the package are available as example data sets.


```{r plotbiomass1 , fig.width=7, fig.height=7, fig.align='center'}
# Load the 1970-1999 version of the North Sea model supplied with the package, run, and generate a plot:
model <- e2e_read("North_Sea", "1970-1999")
results <- e2e_run(model, nyears=3)
e2e_plot_biomass(model, results=results)
```

**Figure 46**. Annual average biomass densities for the final year of a single run of the 1970-1999 North Sea model.

\

```{r plotbiomass2 , fig.width=7, fig.height=7, fig.align='center'}
# For the same model as supplied with the package, plot the example data with credible intervals:
e2e_plot_biomass(model, ci.data=TRUE, use.example=TRUE)
```

**Figure 47**. Annual average biomass densities with credible intervals for the final year of a single run of the 1970-1999 North Sea model.

\

\

# References

Beddington, J.R. (1975). Mutual interference between parasites or predators and its effect on searching efficiency. *Journal of Animal Ecology*, **51**, 331-340.

Bertsimas, D. & Tsitsiklis, J. (1993). Simulated annealing. *Statistical Science*, **8**, 10-15.

Cerny, V. (1985). A thermodynamic approach to the travelling salesman problem: An efficient simulation. *Journal of Optimization Theory and Applications*, **54**, 41-51.

DeAngelis, D.L., Goldstein, R.A.. & O'Neill, R.V. (1975). A model for trophic interaction. *Ecology*, **56**, 881-892.

Heath, M.R., Speirs, D.C. & Steele, J.H. (2014). Understanding patterns and processes in models of trophic cascades. *Ecology Letters*, **17**, 101-114.

Kirkpatrick, S., Gelett, C.D. & Vecchi, M.P. (1983). Optimisation by simulated annealing. *Science*, **220**, 621-630.

Morris, M.D. (1991). Factorial sampling plans for preliminary computational experiments. *Technometrics*, **33**, 161-174.

Pace, M.L., Cole, J.J., Carpenter, S.R. & Kitchell, J.F. (1999) Trophic cascades revealed in diverse ecosystems. *Trends in Ecology and Evolution*, **14**, 483-488.


