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

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.

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.

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

NOAA-EDAB/Rpath documentation built on June 5, 2024, 9:03 p.m.