In zanbi/NOAA-EDAB-Rpath: R implementation of Ecopath with Ecosim

Rpath is an implementation of the ecosystem model Ecopath with Ecosim (EwE; Christensen and Pauly 1992^[Christensen and Pauly. 1992. ECOPATH II - a software for balancing steady-state models and calculating network characteristics. Ecological Modelling 61:169-85], Walters et al. 1997^[Walters et al. 1997. Structuring dynamic models of exploited ecosystems from trophic mass-balance assessments. Reviews of Fish Biology and Fisheries 7:1-34]). This vignette describes some of the basic functionality of the package using a fictional ecosystem, R Ecosystem. Any resemblance to an actual ecosystem is purely coincidental. To see the underlying mathematics please refer to Lucey et al. (in prep^[Lucey et al. in prep. Improving the EBFM toolbox with an alternative open source version of Ecopath with Ecosim]).

knitr::opts_chunk$set(
  comment = '#>',
  collapse = T)
library(Rpath); library(data.table)
````

## Setting up Ecopath

### 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 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.

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.

#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.

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.

#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]

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.

#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)]

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.

REco.params <- rpath.stanzas(REco.params)

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.

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:

REco.params$diet[Group == 'OtherGroundfish', Seabirds := 0.1]

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

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]

knitr::kable(REco.params$diet[, list(Group, Seabirds, Whales)])

Here is the completed model parameter file for R Ecosystem:

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)]

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.

knitr::kable(REco.params$pedigree, caption = 'Pedigree parameters for R Ecosystem')

Running Ecopath

After creating the parameter object, running ecopath in R is relatively straightforward. It is just the function rpath supplied with the parameter object. Additionally, you can supply an ecosystem name for the output.

REco <- rpath(REco.params, eco.name = 'R Ecosystem')
REco

The output object from rpath is an S3 object type called 'Rpath'. Rpath objects are a list of parameters from the mass balance. However, the print function will display the same information as the "Basic Estimates" tab from EwE. You will also notice that the print function will display whether the model is balanced or not. If the model was not balanced, it would list the groups that are not balanced.

You can also display the mortalities associated with each group by supplying the argument morts = T to the print function.

print(REco, morts = T)

Note that if you wish to save the print output you need to use the function write.rpath. This function will also accept the argument 'morts = T'.

The generic function summary will display some summary statistics on the model as well as a list of attributes you can access. To access any of the other attributes simply use the standard list notation.

summary(REco)
REco$TL

One of the advantages of R is its graphical ability. Users can feel free to develop their own graphical routines for the Rpath outputs. However, we have included a basic food web plot. The routine can include fisheries, display group numbers or names, and even highlight a particular group.

webplot(REco)
webplot(REco, labels = T)
webplot(REco, fleets = T, highlight = 'AduRoundfish1')

zanbi/NOAA-EDAB-Rpath documentation built on March 11, 2021, 12:28 a.m.