---
title: "Set up a food web model"
author: "Vignette Author"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Set up a food web model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```


## Setting up a food web model

First, install Rpath as instructed in [Getting Started](). R package `data.table` is also required to run the code in this article.

```{r rpath load}
library(Rpath); library(data.table)
```


### Parameter file generation
Unlike the GUI based EwE software package, Rpath relies on a parameter input file.
This file is actually a list of several different parameter files: model, diet, 
stanzas, and pedigree.  Parameter files can be created outside of R and read in using
the `read.rpath.params` function.  This function will merge several different flat
files into an R object of the list type.  A preferred alternative is to generate 
the list file and populate it completely within R.  The function 
`create.rpath.params` will generate an Rpath.param.  This ensures that all 
of the correct columns are present in the parameter file.

The parameter file contains all of the information you would normally enter in the
input data tabs in EwE.  There are 2 necessary pieces of information to generate
the parameter file: the group names and their corresponding type.  The types are:
living = 0, primary producer = 1, detritus = 2, and fleet = 3.  If your model 
contains multi-stanza groups then you need 2 additional pieces of information: 
stanza group names (include NA for those groups not in a stanza) and the number of
stanzas per stanza group.

```{r groups}
#Groups and types for the R Ecosystem

groups <- c('Seabirds', 'Whales', 'Seals', 'JuvRoundfish1', 'AduRoundfish1', 
            'JuvRoundfish2', 'AduRoundfish2', 'JuvFlatfish1', 'AduFlatfish1',
            'JuvFlatfish2', 'AduFlatfish2', 'OtherGroundfish', 'Foragefish1',
            'Foragefish2', 'OtherForagefish', 'Megabenthos', 'Shellfish',
            'Macrobenthos', 'Zooplankton', 'Phytoplankton', 'Detritus', 
            'Discards', 'Trawlers', 'Midwater', 'Dredgers')

types  <- c(rep(0, 19), 1, rep(2, 2), rep(3, 3))

stgroups <- c(rep(NA, 3), rep('Roundfish1', 2), rep('Roundfish2', 2), 
              rep('Flatfish1', 2), rep('Flatfish2', 2), rep(NA, 14))

REco.params <- create.rpath.params(group = groups, type = types, stgroup = stgroups)
```

REco.params now contains a list of 4 objects: model, diet, stanzas, and pedigree.
The majority of the parameters are populated with NA save those that have logical
default vaules (i.e 0.66667 for VBGF_d).  

### Model parameters
The model parameter list contains the biomass, production to biomass, consumption to 
biomass, etc. parameters as well as the detrital fate parameters and fleet landings and 
discards. 

```{r blank modfile table, echo=FALSE, results='asis'}
knitr::kable(REco.params$model, caption = 'Example of the model list created using the 
             `create.rpath.param` function')
```

Each of the parameter lists are data tables (With the exception of the stanzas list
which is itself a list of an integer and two data tables). Data tables are an 
extension of the classic data frame class.  Advantages of data tables include
simplified indexing which eases the process of populating the parameters.  For example
you can add data to a specific slot or fill an entire column.

```{r How to fill}
#Example of filling specific slots
REco.params$model[Group %in% c('Seals', 'Megabenthos'), EE := 0.8]

#Example of filling an entire column
biomass <- c(0.0149, 0.454, NA, NA, 1.39, NA, 5.553, NA, 5.766, NA,
             0.739, 7.4, 5.1, 4.7, 5.1, NA, 7, 17.4, 23, 10, rep(NA, 5))
REco.params$model[, Biomass := biomass]

```

Note the use of the operator ':=' to assign values.  This is unique to data tables.

```{r Model Table partial, echo = F}
knitr::kable(REco.params$model[, list(Group, Type, Biomass, EE)], 
             caption = 'Example of assigning a specific slot or a whole column')

```

Here are the rest of the columns for the model list.
```{r Model Table sans stanzas}
#Model
biomass <- c(0.0149, 0.454, NA, NA, 1.39, NA, 5.553, NA, 5.766, NA,
             0.739, 7.4, 5.1, 4.7, 5.1, NA, 7, 17.4, 23, 10, rep(NA, 5))

pb <- c(0.098, 0.031, 0.100, 2.026, 0.42, 2.1, 0.425, 1.5, 0.26, 1.1, 0.18, 0.6,
        0.61, 0.65, 1.5, 0.9, 1.3, 7, 39, 240, rep(NA, 5))

qb <- c(76.750, 6.976, 34.455, NA, 2.19, NA, 3.78, NA, 1.44, NA, 1.69,
        1.764, 3.52, 5.65, 3.6, 2.984, rep (NA, 9))

REco.params$model[, Biomass := biomass]
REco.params$model[, PB := pb]
REco.params$model[, QB := qb]

#EE for groups w/o biomass
REco.params$model[Group %in% c('Seals', 'Megabenthos'), EE := 0.8]

#Production to Consumption for those groups without a QB
REco.params$model[Group %in% c('Shellfish', 'Zooplankton'), ProdCons:= 0.25]
REco.params$model[Group == 'Macrobenthos', ProdCons := 0.35]

#Biomass accumulation and unassimilated consumption
REco.params$model[, BioAcc  := c(rep(0, 22), rep(NA, 3))]
REco.params$model[, Unassim := c(rep(0.2, 18), 0.4, rep(0, 3), rep(NA, 3))]

#Detrital Fate
REco.params$model[, Detritus := c(rep(1, 20), rep(0, 5))]
REco.params$model[, Discards := c(rep(0, 22), rep(1, 3))]

#Fisheries
#Landings
trawl  <- c(rep(0, 4), 0.08, 0, 0.32, 0, 0.09, 0, 0.05, 0.2, rep(0, 10), rep(NA, 3))
mid    <- c(rep(0, 12), 0.3, 0.08, 0.02, rep(0, 7), rep(NA, 3))
dredge <- c(rep(0, 15), 0.1, 0.5, rep(0, 5), rep(NA, 3))
REco.params$model[, Trawlers := trawl]
REco.params$model[, Midwater := mid]
REco.params$model[, Dredgers := dredge]

#Discards
trawl.d  <- c(1e-5, 1e-7, 0.001, 0.001, 0.005, 0.001, 0.009, 0.001, 0.04, 0.001,
              0.01, 0.08, 0.001, 0.001, 0.001, rep(0, 7), rep(NA, 3))
mid.d    <- c(rep(0, 2), 0.001, 0.001, 0.01, 0.001, 0.01, rep(0, 4), 0.05, 0.05,
              0.01, 0.01, rep(0, 7), rep(NA, 3))
dredge.d <- c(rep(0, 3), 0.001, 0.05, 0.001, 0.05, 0.001, 0.05, 0.001, 0.01, 0.05,
              rep(0, 3), 0.09, 0.01, 1e-4, rep(0, 4), rep(NA, 3))
REco.params$model[, Trawlers.disc := trawl.d]
REco.params$model[, Midwater.disc := mid.d]
REco.params$model[, Dredgers.disc := dredge.d]
```

```{r Model Table final, echo = F}
knitr::kable(REco.params$model, 
             caption = 'Example of completed model list')

```

###Stanza Parameters
You may have noticed that the biomass and consumption to biomass parameters are
missing from some of the multistanza groups.  Similar to EwE, Rpath calculates those
parameters to ensure that stanza groups support one another (Christensen and Walters
2004^[Christensen and Walters. 2004. Ecopath with Ecosim: methods, capabilities and 
limitations. Ecological Modelling 172:109-139]).  In order to do this, you need to
populate the stanza list.  As mentioned earlier, this is actually a list itself 
containing 3 things: the number of stanza groups, stanza group parameters, and 
individual stanza parameters.  The number of stanzas is automatically populated.  For 
stanza groups you need their von Bertalanffy growth function specialized K and weight at 
50% maturity divided by their weight infinity (relative weight at maturity).  Individual 
stanzas need the first and last month the species is in the stanza, the total mortality 
(Z) on the stanza, and whether or not it is the leading stanza.

```{r Stanza parameters}
#Group parameters
REco.params$stanzas$stgroups[, VBGF_Ksp := c(0.145, 0.295, 0.0761, 0.112)]
REco.params$stanzas$stgroups[, Wmat     := c(0.0769, 0.561, 0.117,  0.321)]

#Individual stanza parameters
REco.params$stanzas$stindiv[, First   := c(rep(c(0, 24), 3), 0, 48)]
REco.params$stanzas$stindiv[, Last    := c(rep(c(23, 400), 3), 47, 400)]
REco.params$stanzas$stindiv[, Z       := c(2.026, 0.42, 2.1, 0.425, 1.5, 
                                           0.26, 1.1, 0.18)]
REco.params$stanzas$stindiv[, Leading := rep(c(F, T), 4)]

```
```{r Stanza Table initial, echo = F}
knitr::kable(REco.params$stanzas$stgroups)
knitr::kable(REco.params$stanzas$stindiv)
```

The final month of the ultimate stanza can be set to any value.  The function 
`rpath.stanzas` will calculate the final month as the point where the species reaches 
90% Winf.  The function `rpath.stanzas` will also add data tables containing the weight, number, and consumption at age for each stanza group.  

```{r rpath.stanzas}
REco.params <- rpath.stanzas(REco.params)
```
```{r Stanza Table final, echo = F}
knitr::kable(REco.params$stanzas$stanzas, caption = 'Completed stanzas table')
knitr::kable(head(REco.params$stanzas$StGroup[[1]]), 
             caption = 'Example of the StGroup data table')
```

Output from the `rpath.stanzas` function can be plotted using the `stanzaplot`
function.
```{r stanzaplot, fig.align = 'center', fig.height = 5, fig.width = 9}
stanzaplot(REco.params, StanzaGroup = 1)
```
Note: If you do not have multistanza groups in your model, you do not have to run
`rpath.stanzas`.

###Diet Parameters
The data entered in the diet list is the same as the data entered in the 
diet composition tab in EwE. Just as within EwE, the columns 
represent the predators while the rows represent the prey.  Individual diet components 
can be adjusted by specifying the prey in the 'Group' variable and assigning a value to 
the predator.  For example, if you wanted to assign 10% of the seabird diet as 'Other Groundfish' you could do it like this:

```{r how to fill diet 1}
REco.params$diet[Group == 'OtherGroundfish', Seabirds := 0.1]
```

You can also assign the entire diet composition for a predator:

```{r how to fill diet 2}
whale.diet <- c(rep(NA, 3), 0.01, NA, 0.01, NA, 0.01, NA, 0.01, rep(NA, 4), 0.1,
                rep(NA, 3), 0.86, rep(NA, 3), NA)
REco.params$diet[, Whales := whale.diet]
```
```{r Dietfile table partial, echo = F}
knitr::kable(REco.params$diet[, list(Group, Seabirds, Whales)])
```

Here is the completed model parameter file for R Ecosystem:

```{r diet fill}
REco.params$diet[, Seabirds        := c(rep(NA, 11), 0.1, 0.25, 0.2, 0.15, 
                                        rep(NA, 6), 0.3, NA)]
REco.params$diet[, Whales          := c(rep(NA, 3), 0.01, NA, 0.01, NA, 0.01, 
                                        NA, 0.01, rep(NA, 4), 0.1, rep(NA, 3), 
                                        0.86, rep(NA, 3), NA)]
REco.params$diet[, Seals           := c(rep(NA, 3), 0.05, 0.1, 0.05, 0.2, 0.005, 
                                        0.05, 0.005, 0.01, 0.24, rep(0.05, 4), 
                                        0.09, rep(NA, 5), NA)]
REco.params$diet[, JuvRoundfish1   := c(rep(NA, 3), rep(c(1e-4, NA), 4), 1e-3, 
                                        rep(NA, 2), 0.05, 1e-4, NA, .02, 0.7785, 
                                        0.1, 0.05, NA, NA)]
REco.params$diet[, AduRoundfish1   := c(rep(NA, 5), 1e-3, 0.01, 1e-3, 0.05, 1e-3, 
                                        0.01, 0.29, 0.1, 0.1, 0.347, 0.03, NA, 
                                        0.05, 0.01, rep(NA, 3), NA)]
REco.params$diet[, JuvRoundfish2   := c(rep(NA, 3), rep(c(1e-4, NA), 4), 1e-3, 
                                        rep(NA, 2), 0.05, 1e-4, NA, .02, 0.7785, 
                                        0.1, .05, NA, NA)]
REco.params$diet[, AduRoundfish2   := c(rep(NA, 3), 1e-4, NA, 1e-4, NA, rep(1e-4, 4), 
                                        0.1, rep(0.05, 3), 0.2684, 0.01, 0.37, 0.001, 
                                        NA, 0.1, NA, NA)]
REco.params$diet[, JuvFlatfish1    := c(rep(NA, 3), rep(c(1e-4, NA), 4), rep(NA, 3), 
                                        rep(1e-4, 2), NA, 0.416, 0.4334, 0.1, 0.05, 
                                        NA, NA)]
REco.params$diet[, AduFlatfish1    := c(rep(NA, 7), rep(1e-4, 5), rep(NA, 2), 0.001, 
                                        0.05, 0.001, 0.6, 0.2475, NA, 0.1, NA, NA)]
REco.params$diet[, JuvFlatfish2    := c(rep(NA, 3), rep(c(1e-4, NA), 4), rep(NA, 3),
                                        rep(1e-4, 2), NA, 0.416, 0.4334, 0.1, 0.05, 
                                        NA, NA)]
REco.params$diet[, AduFlatfish2    := c(rep(NA, 7), 1e-4, NA, 1e-4, rep(NA, 4), 
                                        rep(1e-4, 3), 0.44, 0.3895, NA, 0.17, NA, NA)]
REco.params$diet[, OtherGroundfish := c(rep(NA, 3), rep(1e-4, 8), 0.05, 0.08, 0.0992, 
                                        0.3, 0.15, 0.01, 0.3, 0.01, rep(NA, 3), NA)]
REco.params$diet[, Foragefish1     := c(rep(NA, 3), rep(c(1e-4, NA), 4), rep(NA, 7), 
                                        0.8196, 0.06, 0.12, NA, NA)]
REco.params$diet[, Foragefish2     := c(rep(NA, 3), rep(c(1e-4, NA), 4), rep(NA, 7), 
                                        0.8196, 0.06, 0.12, NA, NA)]
REco.params$diet[, OtherForagefish := c(rep(NA, 3), rep(c(1e-4, NA), 4), rep(NA, 7), 
                                        0.8196, 0.06, 0.12, NA, NA)]
REco.params$diet[, Megabenthos     := c(rep(NA, 15), 0.1, 0.03, 0.55, rep(NA, 2), 0.32,
                                        NA, NA)]
REco.params$diet[, Shellfish       := c(rep(NA, 18), 0.3, 0.5, 0.2, NA, NA)]
REco.params$diet[, Macrobenthos    := c(rep(NA, 16), 0.01, rep(0.2, 2), NA, 0.59, NA, NA)]
REco.params$diet[, Zooplankton     := c(rep(NA, 18), 0.2, 0.6, 0.2, NA, NA)]
```
```{r Dietfile Table, echo = F}
knitr::kable(REco.params$diet, caption = 'Diet parameters for R Ecosystem')
```

###Pedigree parameters
Rpath does not currently use pedigrees however, future Rpath extensions will use
them.  Therefore we include them in the current parameter object. The
default values are 1 (low confidence).  These defaults are not changed for R
Ecosystem but can obviously be changed in a similar manner to the other parameter
files.
```{r pedigree table, echo = F}
knitr::kable(REco.params$pedigree, caption = 'Pedigree parameters for R Ecosystem')
```